1. Executive Summary

Comprehensive regional development rooted in multi-nodal growth engines has been the new norm in the context of American urban planning. Cross-boundary collaboration brings benefits towards resource distribution, a balanced growth and overall sustainable development. Conjuncting Pennsylvania, New Jersey, and Delaware and in proximity to New York City and Washington D.C., the Delaware Valley Region (DVR) serves as an example of regional collaboration in growth management. With the forecasts of next ten-to-thirty years development, the region witnesses the potential of increased population joining the workforce, along with the expansion and enhancement of regional infrastructure networks. This report will summarize our recent study predicting future development demand based on the demand-driven factor of population change and the supply-driven factor of infrastructure network enhancement. It is predicted that development would intensify in Chester and Gloucester County as the northern parts of Montgomery, Bucks and Burlington County. The development would be further intensified along the proposed highway, concentrating in north Montgomery and Burlington County. The findings provide insights for planners to formulate better policies to manage urban sprawl, leverage smart growth and balance regional development.

1.2 Setup map themes

library(tidyverse)
## ── Attaching packages ───────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.0     ✓ purrr   0.3.4
## ✓ tibble  3.0.1     ✓ dplyr   0.8.5
## ✓ tidyr   1.0.3     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## ── Conflicts ──────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(sf)
## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
library(raster)
## Loading required package: sp
## 
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
## 
##     select
## The following object is masked from 'package:tidyr':
## 
##     extract
library(knitr)
library(kableExtra)
## 
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows
library(tidycensus)
library(tigris)
## To enable 
## caching of data, set `options(tigris_use_cache = TRUE)` in your R script or .Rprofile.
## 
## Attaching package: 'tigris'
## The following object is masked from 'package:graphics':
## 
##     plot
## The following object is masked from 'package:base':
## 
##     plot
library(FNN)
library(QuantPsyc)
## Loading required package: boot
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following objects are masked from 'package:raster':
## 
##     area, select
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## Attaching package: 'QuantPsyc'
## The following object is masked from 'package:base':
## 
##     norm
library(caret)
## Loading required package: lattice
## 
## Attaching package: 'lattice'
## The following object is masked from 'package:boot':
## 
##     melanoma
## 
## Attaching package: 'caret'
## The following object is masked from 'package:purrr':
## 
##     lift
library(yardstick)
## For binary classification, the first factor level is assumed to be the event.
## Set the global option `yardstick.event_first` to `FALSE` to change this.
## 
## Attaching package: 'yardstick'
## The following objects are masked from 'package:caret':
## 
##     precision, recall, sensitivity, specificity
## The following object is masked from 'package:readr':
## 
##     spec
library(pscl)
## Classes and Methods for R developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University
## Simon Jackman
## hurdle and zeroinfl functions by Achim Zeileis
library(plotROC) 
library(ggrepel)
library(pROC)
## Type 'citation("pROC")' for a citation.
## 
## Attaching package: 'pROC'
## The following object is masked from 'package:plotROC':
## 
##     ggroc
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
library(grid)
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
library(viridis)
## Loading required package: viridisLite
# set up map themes
mapTheme <- function(base_size = 12) {
  theme(
    text = element_text( color = "black"),
    plot.title = element_text(size = 14,colour = "black"),
    plot.subtitle=element_text(face="italic"),
    plot.caption=element_text(hjust=0),
    axis.ticks = element_blank(),
    panel.background = element_blank(),axis.title = element_blank(),
    axis.text = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(colour = "black", fill=NA, size=2)
  )
}

plotTheme <- function(base_size = 12) {
  theme(
    text = element_text( color = "black"),
    plot.title = element_text(size = 14,colour = "black"),
    plot.subtitle = element_text(face="italic"),
    plot.caption = element_text(hjust=0),
    axis.ticks = element_blank(),
    panel.background = element_blank(),
    panel.grid.major = element_line("grey80", size = 0.1),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(colour = "black", fill=NA, size=2),
    strip.background = element_rect(fill = "grey80", color = "white"),
    strip.text = element_text(size=12),
    axis.title = element_text(size=12),
    axis.text = element_text(size=10),
    plot.background = element_blank(),
    legend.background = element_blank(),
    legend.title = element_text(colour = "black", face = "italic"),
    legend.text = element_text(colour = "black", face = "italic"),
    strip.text.x = element_text(size = 14)
  )
}

palette2 <- c("#41b6c4","#253494")
palette4 <- c("#a1dab4","#41b6c4","#2c7fb8","#253494")
palette5 <- c("#ffffcc","#a1dab4","#41b6c4","#2c7fb8","#253494")
palette10 <- c("#f7fcf0","#e0f3db","#ccebc5","#a8ddb5","#7bccc4",
               "#4eb3d3","#2b8cbe","#0868ac","#084081","#f7fcf0")

#this function converts a column in to quintiles. It is used for mapping.
quintileBreaks <- function(df,variable) {
  as.character(quantile(df[[variable]],
                        c(.01,.2,.4,.6,.8),na.rm=T))
}

#This function can be used to convert a polygon sf to centroids xy coords.
xyC <- function(aPolygonSF) {
  as.data.frame(
    cbind(x=st_coordinates(st_centroid(aPolygonSF))[,1],
          y=st_coordinates(st_centroid(aPolygonSF))[,2]))
} 

#this function convert a raster to a data frame so it can be plotted in ggplot
rast <- function(inRaster) {
  data.frame(
    xyFromCell(inRaster, 1:ncell(inRaster)), 
    value = getValues(inRaster)) }
  1. Data Wrangling and feature engineering 2.2 Land cover change data
DVR <- 
  st_read("DVR_boundaryproj.shp") %>%
  st_transform(102729)
## Reading layer `DVR_boundaryproj' from data source `/Users/haitianwang/Documents/675/Final/DVR_boundaryproj.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 2422314 ymin: 78345.53 xmax: 2915049 ymax: 474801.3
## CRS:            2272
lc_change = raster("lc_chang_1kproj.tif")

reclassMatrix <- 
  matrix(c(
    0,0,0,
    1,1,1,
    1,Inf,0),
    ncol=3, byrow=T)

reclassMatrix
##      [,1] [,2] [,3]
## [1,]    0    0    0
## [2,]    1    1    1
## [3,]    1  Inf    0
lc_change2 <- 
  reclassify(lc_change,reclassMatrix)

lc_change2[lc_change2 < 1] <- NA

names(lc_change2) <- "lc_change"

ggplot() +
  geom_sf(data=DVR) +
  geom_raster(data=rast(lc_change2) %>% na.omit, 
              aes(x,y,fill=as.factor(value))) +
  scale_fill_viridis(discrete=TRUE, name ="Land Cover\nChange") + 
  labs(title="Development land use change") +
  mapTheme()
## Warning: Raster pixels are placed at uneven horizontal intervals and will be
## shifted. Consider using geom_tile() instead.
## Warning: Raster pixels are placed at uneven vertical intervals and will be
## shifted. Consider using geom_tile() instead.

ggplot() +
  geom_sf(data=DVR) +
  geom_raster(data=rast(lc_change) %>% na.omit, 
              aes(x,y,fill=as.factor(value))) +
  scale_fill_viridis(discrete=TRUE, name ="Land Cover\nChange") + 
  labs(title="Development land use change") +
  mapTheme()
## Warning: Raster pixels are placed at uneven horizontal intervals and will be
## shifted. Consider using geom_tile() instead.

Create fishnet

DVR_fishnet <- 
  st_make_grid(DVR, 4000) %>%
  st_sf()

DVR_fishnet <-
  DVR_fishnet[DVR,]

ggplot() +
  geom_sf(data=DVR_fishnet) +
  labs(title="Fishnet, 3000 foot resolution") +
  mapTheme()

changePoints <-
  rasterToPoints(lc_change2) %>%
  as.data.frame() %>%
  st_as_sf(coords = c("x", "y"), crs = st_crs(DVR_fishnet))




ggplot() +
  geom_sf(data=changePoints) +
  labs(title="Fishnet, 3000 foot resolution") +
  mapTheme()

st_crs(changePoints) == st_crs(DVR_fishnet)
## [1] TRUE
ggplot() +
  geom_sf(data=DVR_fishnet, fill="black", colour=NA) +
  geom_sf(data=filter(changePoints, lc_change == 1), colour = "white")

head(DVR_fishnet)
## Simple feature collection with 6 features and 0 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 2738314 ymin: 78345.53 xmax: 2754314 ymax: 86345.53
## CRS:            2272
##                         geometry
## 1 POLYGON ((2742314 78345.53,...
## 2 POLYGON ((2746314 78345.53,...
## 3 POLYGON ((2750314 78345.53,...
## 4 POLYGON ((2738314 82345.53,...
## 5 POLYGON ((2742314 82345.53,...
## 6 POLYGON ((2746314 82345.53,...
fishnet <- 
  aggregate(changePoints, DVR_fishnet, sum) %>%
  mutate(lc_change = ifelse(is.na(lc_change),0,1),
         lc_change = as.factor(lc_change))

plot(fishnet)

str(fishnet)
## Classes 'sf' and 'data.frame':   6934 obs. of  2 variables:
##  $ lc_change: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
##  $ geometry :sfc_POLYGON of length 6934; first list element: List of 1
##   ..$ : num [1:5, 1:2] 2742314 2746314 2746314 2742314 2742314 ...
##   ..- attr(*, "class")= chr [1:3] "XY" "POLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA
##   ..- attr(*, "names")= chr "lc_change"
na.omit(fishnet, lc_change) %>% 
ggplot() +
  geom_sf(data=DVR) +
  geom_point(data=fishnet, 
             aes(x=xyC(fishnet)$x, y=xyC(fishnet)$y, colour=lc_change)) +
    scale_colour_manual(values = palette2,
                      labels=c("No Change","New Development"),
                      name = "") +
  labs(title = "Land cover development change", subtitle = "As fishnet centroids") +
  mapTheme()
## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

#rlang::last_error()

2.3 Land Cover in 2001

lc_2001 <- raster("lc_2001_1kproj.tif")

ggplot() +
  geom_sf(data=DVR) +
  geom_raster(data=rast(lc_2001) %>% na.omit %>% filter(value > 0), 
              aes(x,y,fill=as.factor(value))) +
  scale_fill_viridis(discrete=TRUE, name ="") +
  labs(title = "Land Cover, 2001") +
  mapTheme()
## Warning: Raster pixels are placed at uneven horizontal intervals and will be
## shifted. Consider using geom_tile() instead.

developed <- lc_2001 == 21 | lc_2001 == 22 | lc_2001 == 23 | lc_2001 == 24
forest <- lc_2001 == 41 | lc_2001 == 42 | lc_2001 == 43 
farm <- lc_2001 == 81 | lc_2001 == 82 
wetlands <- lc_2001 == 90 | lc_2001 == 95 
otherUndeveloped <- lc_2001 == 52 | lc_2001 == 71 | lc_2001 == 31 
water <- lc_2001 == 11

names(developed) <- "developed"
names(forest) <- "forest"
names(farm) <- "farm"
names(wetlands) <- "wetlands"
names(otherUndeveloped) <- "otherUndeveloped"
names(water) <- "water"



aggregateRaster <- function(inputRasterList, theFishnet) {
  #create an empty fishnet with the same dimensions as the input fishnet
  theseFishnets <- theFishnet %>% dplyr::select()
  #for each raster in the raster list
  for (i in inputRasterList) {
    #create a variable name corresponding to the ith raster
    varName <- names(i)
    #convert raster to points as an sf
    thesePoints <-
      rasterToPoints(i) %>%
      as.data.frame() %>%
      st_as_sf(coords = c("x", "y"), crs = st_crs(theFishnet)) %>%
      filter(.[[1]] == 1)
    #aggregate to the fishnet
    thisFishnet <-
      aggregate(thesePoints, theFishnet, length) %>%
      mutate(!!varName := ifelse(is.na(.[[1]]),0,1))
    #add to the larger fishnet
    theseFishnets <- cbind(theseFishnets,thisFishnet)
  }
  #output all aggregates as one large fishnet
  return(theseFishnets)
}





theRasterList <- c(developed,forest,farm,wetlands,otherUndeveloped,water)

aggregatedRasters <-
  aggregateRaster(theRasterList, DVR_fishnet) %>%
  dplyr::select(developed,forest,farm,wetlands,otherUndeveloped,water) %>%
  mutate_if(is.numeric,as.factor)


aggregatedRasters %>%
  gather(var,value,developed:water) %>%
  st_cast("POLYGON") %>%    #just to make sure no weird geometries slipped in
  mutate(X = xyC(.)$x,
         Y = xyC(.)$y) %>%
  ggplot() +
  geom_sf(data=DVR) +
  geom_point(aes(X,Y, colour=as.factor(value))) +
  facet_wrap(~var) +
  scale_colour_manual(values = palette2, na.translate = F,
                      labels=c("Other","Land Cover"),
                      name = "") +
  labs(title = "Land cover types, 2001",
       subtitle = "As fishnet centroids") +
  mapTheme()
## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

# 2.4 Census Data

census_api_key("95e1ba2797ea00d4a00a86bf8364df089c47a047", install = TRUE , overwrite=TRUE)
## Your original .Renviron will be backed up and stored in your R HOME directory if needed.
## Your API key has been stored in your .Renviron and can be accessed by Sys.getenv("CENSUS_API_KEY"). 
## To use now, restart R or run `readRenviron("~/.Renviron")`
## [1] "95e1ba2797ea00d4a00a86bf8364df089c47a047"
PAPop00 <- 
  get_decennial(geography = "tract", variables = "P001001", year = 2000,
                state =42, geometry = TRUE, 
                county=c("Bucks","Montgomery",
                         "Chester","Philadelphia","Delaware"))
## Getting data from the 2000 decennial Census
## Downloading feature geometry from the Census website.  To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=                                                                     |   1%
  |                                                                            
  |==                                                                    |   2%
  |                                                                            
  |==                                                                    |   3%
  |                                                                            
  |===                                                                   |   4%
  |                                                                            
  |===                                                                   |   5%
  |                                                                            
  |====                                                                  |   6%
  |                                                                            
  |=====                                                                 |   7%
  |                                                                            
  |=====                                                                 |   8%
  |                                                                            
  |======                                                                |   9%
  |                                                                            
  |=======                                                               |   9%
  |                                                                            
  |=======                                                               |  11%
  |                                                                            
  |========                                                              |  12%
  |                                                                            
  |=========                                                             |  12%
  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |===========                                                           |  15%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |=============                                                         |  18%
  |                                                                            
  |=============                                                         |  19%
  |                                                                            
  |==============                                                        |  20%
  |                                                                            
  |===============                                                       |  21%
  |                                                                            
  |===============                                                       |  22%
  |                                                                            
  |================                                                      |  23%
  |                                                                            
  |=================                                                     |  24%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |==================                                                    |  26%
  |                                                                            
  |===================                                                   |  27%
  |                                                                            
  |====================                                                  |  28%
  |                                                                            
  |====================                                                  |  29%
  |                                                                            
  |=====================                                                 |  30%
  |                                                                            
  |======================                                                |  31%
  |                                                                            
  |======================                                                |  32%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |========================                                              |  35%
  |                                                                            
  |=========================                                             |  36%
  |                                                                            
  |==========================                                            |  38%
  |                                                                            
  |===========================                                           |  38%
  |                                                                            
  |============================                                          |  39%
  |                                                                            
  |============================                                          |  41%
  |                                                                            
  |=============================                                         |  41%
  |                                                                            
  |==============================                                        |  42%
  |                                                                            
  |==============================                                        |  43%
  |                                                                            
  |===============================                                       |  44%
  |                                                                            
  |================================                                      |  45%
  |                                                                            
  |================================                                      |  46%
  |                                                                            
  |=================================                                     |  47%
  |                                                                            
  |=================================                                     |  48%
  |                                                                            
  |==================================                                    |  49%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |===================================                                   |  51%
  |                                                                            
  |====================================                                  |  52%
  |                                                                            
  |=====================================                                 |  52%
  |                                                                            
  |======================================                                |  54%
  |                                                                            
  |======================================                                |  55%
  |                                                                            
  |=======================================                               |  55%
  |                                                                            
  |========================================                              |  57%
  |                                                                            
  |=========================================                             |  58%
  |                                                                            
  |==========================================                            |  60%
  |                                                                            
  |===========================================                           |  61%
  |                                                                            
  |===========================================                           |  62%
  |                                                                            
  |============================================                          |  63%
  |                                                                            
  |=============================================                         |  64%
  |                                                                            
  |==============================================                        |  65%
  |                                                                            
  |==============================================                        |  66%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |================================================                      |  68%
  |                                                                            
  |================================================                      |  69%
  |                                                                            
  |=================================================                     |  70%
  |                                                                            
  |==================================================                    |  71%
  |                                                                            
  |==================================================                    |  72%
  |                                                                            
  |===================================================                   |  73%
  |                                                                            
  |====================================================                  |  74%
  |                                                                            
  |====================================================                  |  75%
  |                                                                            
  |=====================================================                 |  76%
  |                                                                            
  |======================================================                |  78%
  |                                                                            
  |=======================================================               |  79%
  |                                                                            
  |========================================================              |  79%
  |                                                                            
  |========================================================              |  81%
  |                                                                            
  |=========================================================             |  81%
  |                                                                            
  |==========================================================            |  82%
  |                                                                            
  |==========================================================            |  84%
  |                                                                            
  |===========================================================           |  84%
  |                                                                            
  |============================================================          |  85%
  |                                                                            
  |============================================================          |  86%
  |                                                                            
  |=============================================================         |  87%
  |                                                                            
  |==============================================================        |  88%
  |                                                                            
  |==============================================================        |  89%
  |                                                                            
  |===============================================================       |  90%
  |                                                                            
  |===============================================================       |  91%
  |                                                                            
  |================================================================      |  92%
  |                                                                            
  |=================================================================     |  93%
  |                                                                            
  |==================================================================    |  94%
  |                                                                            
  |==================================================================    |  95%
  |                                                                            
  |===================================================================   |  95%
  |                                                                            
  |====================================================================  |  97%
  |                                                                            
  |===================================================================== |  98%
  |                                                                            
  |======================================================================| 100%
NJPop00 <- 
  get_decennial(geography = "tract", variables = "P001001", year = 2000,
                state =34, geometry = TRUE, 
                county=c("Mercer","Camden","Gloucester","Burlington"))  
## Getting data from the 2000 decennial Census
## Downloading feature geometry from the Census website.  To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |==                                                                    |   3%
  |                                                                            
  |=====                                                                 |   6%
  |                                                                            
  |=====                                                                 |   7%
  |                                                                            
  |=======                                                               |  10%
  |                                                                            
  |=========                                                             |  13%
  |                                                                            
  |===========                                                           |  16%
  |                                                                            
  |==============                                                        |  19%
  |                                                                            
  |===============                                                       |  21%
  |                                                                            
  |=================                                                     |  24%
  |                                                                            
  |==================                                                    |  26%
  |                                                                            
  |====================                                                  |  29%
  |                                                                            
  |=======================                                               |  32%
  |                                                                            
  |========================                                              |  34%
  |                                                                            
  |==========================                                            |  37%
  |                                                                            
  |===========================                                           |  39%
  |                                                                            
  |=============================                                         |  42%
  |                                                                            
  |================================                                      |  45%
  |                                                                            
  |=================================                                     |  47%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |====================================                                  |  52%
  |                                                                            
  |=======================================                               |  55%
  |                                                                            
  |=========================================                             |  58%
  |                                                                            
  |==========================================                            |  60%
  |                                                                            
  |============================================                          |  63%
  |                                                                            
  |=============================================                         |  65%
  |                                                                            
  |================================================                      |  68%
  |                                                                            
  |==================================================                    |  71%
  |                                                                            
  |===================================================                   |  73%
  |                                                                            
  |=====================================================                 |  76%
  |                                                                            
  |======================================================                |  78%
  |                                                                            
  |=========================================================             |  81%
  |                                                                            
  |===========================================================           |  84%
  |                                                                            
  |============================================================          |  86%
  |                                                                            
  |==============================================================        |  89%
  |                                                                            
  |===============================================================       |  91%
  |                                                                            
  |==================================================================    |  94%
  |                                                                            
  |====================================================================  |  97%
  |                                                                            
  |===================================================================== |  99%
  |                                                                            
  |======================================================================| 100%
DVRPop00 <- rbind(PAPop00, NJPop00)%>%
  rename(pop_2000 = value) %>%
  st_transform(st_crs(DVR_fishnet)) 

 
PAPop10 <- 
  get_decennial(geography = "tract", variables = "P001001", year = 2010,
                state =42, geometry = TRUE, 
                county=c("Bucks","Montgomery",
                         "Chester","Philadelphia","Delaware"))
## Getting data from the 2010 decennial Census
## Downloading feature geometry from the Census website.  To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=                                                                     |   1%
  |                                                                            
  |==                                                                    |   2%
  |                                                                            
  |==                                                                    |   3%
  |                                                                            
  |===                                                                   |   4%
  |                                                                            
  |===                                                                   |   5%
  |                                                                            
  |====                                                                  |   6%
  |                                                                            
  |=====                                                                 |   7%
  |                                                                            
  |=====                                                                 |   8%
  |                                                                            
  |======                                                                |   9%
  |                                                                            
  |=======                                                               |  10%
  |                                                                            
  |========                                                              |  11%
  |                                                                            
  |========                                                              |  12%
  |                                                                            
  |=========                                                             |  13%
  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |==========                                                            |  15%
  |                                                                            
  |===========                                                           |  16%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |============                                                          |  18%
  |                                                                            
  |=============                                                         |  19%
  |                                                                            
  |==============                                                        |  19%
  |                                                                            
  |==============                                                        |  21%
  |                                                                            
  |===============                                                       |  22%
  |                                                                            
  |================                                                      |  22%
  |                                                                            
  |=================                                                     |  24%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |===================                                                   |  27%
  |                                                                            
  |====================                                                  |  29%
  |                                                                            
  |=====================                                                 |  30%
  |                                                                            
  |======================                                                |  32%
  |                                                                            
  |=======================                                               |  32%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |========================                                              |  34%
  |                                                                            
  |=========================                                             |  35%
  |                                                                            
  |==========================                                            |  36%
  |                                                                            
  |==========================                                            |  37%
  |                                                                            
  |===========================                                           |  38%
  |                                                                            
  |===========================                                           |  39%
  |                                                                            
  |============================                                          |  40%
  |                                                                            
  |=============================                                         |  41%
  |                                                                            
  |=============================                                         |  42%
  |                                                                            
  |==============================                                        |  43%
  |                                                                            
  |===============================                                       |  44%
  |                                                                            
  |===============================                                       |  45%
  |                                                                            
  |================================                                      |  46%
  |                                                                            
  |=================================                                     |  47%
  |                                                                            
  |==================================                                    |  48%
  |                                                                            
  |==================================                                    |  49%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |====================================                                  |  51%
  |                                                                            
  |====================================                                  |  52%
  |                                                                            
  |=====================================                                 |  53%
  |                                                                            
  |======================================                                |  55%
  |                                                                            
  |=======================================                               |  56%
  |                                                                            
  |========================================                              |  56%
  |                                                                            
  |========================================                              |  58%
  |                                                                            
  |=========================================                             |  58%
  |                                                                            
  |==========================================                            |  60%
  |                                                                            
  |===========================================                           |  61%
  |                                                                            
  |============================================                          |  63%
  |                                                                            
  |=============================================                         |  64%
  |                                                                            
  |==============================================                        |  66%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |================================================                      |  68%
  |                                                                            
  |================================================                      |  69%
  |                                                                            
  |=================================================                     |  70%
  |                                                                            
  |==================================================                    |  71%
  |                                                                            
  |===================================================                   |  72%
  |                                                                            
  |===================================================                   |  73%
  |                                                                            
  |====================================================                  |  74%
  |                                                                            
  |=====================================================                 |  75%
  |                                                                            
  |=====================================================                 |  76%
  |                                                                            
  |======================================================                |  77%
  |                                                                            
  |======================================================                |  78%
  |                                                                            
  |=======================================================               |  79%
  |                                                                            
  |========================================================              |  80%
  |                                                                            
  |=========================================================             |  81%
  |                                                                            
  |=========================================================             |  82%
  |                                                                            
  |==========================================================            |  83%
  |                                                                            
  |===========================================================           |  84%
  |                                                                            
  |============================================================          |  85%
  |                                                                            
  |============================================================          |  86%
  |                                                                            
  |=============================================================         |  87%
  |                                                                            
  |=============================================================         |  88%
  |                                                                            
  |==============================================================        |  89%
  |                                                                            
  |===============================================================       |  90%
  |                                                                            
  |===============================================================       |  91%
  |                                                                            
  |================================================================      |  92%
  |                                                                            
  |=================================================================     |  92%
  |                                                                            
  |==================================================================    |  94%
  |                                                                            
  |==================================================================    |  95%
  |                                                                            
  |===================================================================   |  95%
  |                                                                            
  |====================================================================  |  97%
  |                                                                            
  |===================================================================== |  98%
  |                                                                            
  |======================================================================| 100%
NJPop10 <- 
  get_decennial(geography = "tract", variables = "P001001", year = 2010,
                state =34, geometry = TRUE, 
                county=c("Mercer","Camden","Gloucester","Burlington"))  
## Getting data from the 2010 decennial Census
## Downloading feature geometry from the Census website.  To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |==                                                                    |   2%
  |                                                                            
  |===                                                                   |   5%
  |                                                                            
  |====                                                                  |   5%
  |                                                                            
  |======                                                                |   8%
  |                                                                            
  |=======                                                               |  10%
  |                                                                            
  |=========                                                             |  12%
  |                                                                            
  |==========                                                            |  15%
  |                                                                            
  |===========                                                           |  16%
  |                                                                            
  |=============                                                         |  18%
  |                                                                            
  |==============                                                        |  20%
  |                                                                            
  |===============                                                       |  22%
  |                                                                            
  |=================                                                     |  24%
  |                                                                            
  |==================                                                    |  26%
  |                                                                            
  |====================                                                  |  28%
  |                                                                            
  |=====================                                                 |  29%
  |                                                                            
  |======================                                                |  32%
  |                                                                            
  |========================                                              |  34%
  |                                                                            
  |=========================                                             |  35%
  |                                                                            
  |===========================                                           |  38%
  |                                                                            
  |===========================                                           |  39%
  |                                                                            
  |=============================                                         |  42%
  |                                                                            
  |===============================                                       |  44%
  |                                                                            
  |================================                                      |  45%
  |                                                                            
  |=================================                                     |  48%
  |                                                                            
  |==================================                                    |  49%
  |                                                                            
  |====================================                                  |  51%
  |                                                                            
  |======================================                                |  54%
  |                                                                            
  |=======================================                               |  55%
  |                                                                            
  |========================================                              |  58%
  |                                                                            
  |=========================================                             |  59%
  |                                                                            
  |===========================================                           |  61%
  |                                                                            
  |=============================================                         |  64%
  |                                                                            
  |=============================================                         |  65%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |================================================                      |  69%
  |                                                                            
  |==================================================                    |  71%
  |                                                                            
  |===================================================                   |  74%
  |                                                                            
  |====================================================                  |  75%
  |                                                                            
  |======================================================                |  77%
  |                                                                            
  |=======================================================               |  78%
  |                                                                            
  |=========================================================             |  81%
  |                                                                            
  |==========================================================            |  83%
  |                                                                            
  |===========================================================           |  85%
  |                                                                            
  |=============================================================         |  87%
  |                                                                            
  |==============================================================        |  88%
  |                                                                            
  |================================================================      |  91%
  |                                                                            
  |=================================================================     |  93%
  |                                                                            
  |==================================================================    |  94%
  |                                                                            
  |====================================================================  |  97%
  |                                                                            
  |===================================================================== |  98%
  |                                                                            
  |======================================================================| 100%
DVRPop10 <- rbind(PAPop10, NJPop10)%>%
  rename(pop_2010 = value) %>%
  st_transform(st_crs(DVR_fishnet)) %>%
  st_buffer(-1)


grid.arrange(
  ggplot() +
    geom_sf(data = DVRPop00, aes(fill=factor(ntile(pop_2000,5))), colour=NA) +
    scale_fill_manual(values = palette5,
                      labels=quintileBreaks(DVRPop00,"pop_2000"),
                      name="Quintile\nBreaks") +
    labs(title="Population, DVR: 2000") +
    mapTheme(),
  
  ggplot() +
    geom_sf(data = DVRPop10, aes(fill=factor(ntile(pop_2010,5))), colour=NA) +
    scale_fill_manual(values = palette5,
                      labels=quintileBreaks(DVRPop10,"pop_2010"),
                      name="Quintile\nBreaks") +
    labs(title="Population, DVR: 2010") +
    mapTheme(), ncol=2)

DVR_fishnet <-
  DVR_fishnet %>%
  rownames_to_column("fishnetID") %>% 
  mutate(fishnetID = as.numeric(fishnetID)) %>%
  dplyr::select(fishnetID)

fishnetPopulation00 <-
  st_interpolate_aw(DVRPop00["pop_2000"],DVR_fishnet, extensive=TRUE) %>%
  as.data.frame(.) %>%
  left_join(DVR_fishnet, ., by=c("fishnetID"='Group.1')) %>% 
  mutate(pop_2000 = replace_na(pop_2000,0)) %>%
  dplyr::select(pop_2000)
## Warning in st_interpolate_aw.sf(DVRPop00["pop_2000"], DVR_fishnet, extensive =
## TRUE): st_interpolate_aw assumes attributes are constant over areas of x
fishnetPopulation10 <-
  st_interpolate_aw(DVRPop10["pop_2010"],DVR_fishnet, extensive=TRUE) %>%
  as.data.frame(.) %>%
  left_join(DVR_fishnet, ., by=c("fishnetID"='Group.1')) %>% 
  mutate(pop_2010 = replace_na(pop_2010,0)) %>%
  dplyr::select(pop_2010)
## Warning in st_interpolate_aw.sf(DVRPop10["pop_2010"], DVR_fishnet, extensive =
## TRUE): st_interpolate_aw assumes attributes are constant over areas of x
fishnetPopulation <- 
  cbind(fishnetPopulation00,fishnetPopulation10) %>%
  dplyr::select(pop_2000,pop_2010) %>%
  mutate(pop_Change = pop_2010 - pop_2000)
grid.arrange(
  ggplot() +
    geom_sf(data=DVRPop10, aes(fill=factor(ntile(pop_2010,5))),colour=NA) +
    scale_fill_manual(values = palette5,
                      labels=substr(quintileBreaks(DVRPop10,"pop_2010"),1,4),
                      name="Quintile\nBreaks") +
    labs(title="Population, DVR: 2010",
         subtitle="Represented as tracts; Boundaries omitted") +
    mapTheme(),
  
  ggplot() +
    geom_sf(data=fishnetPopulation, aes(fill=factor(ntile(pop_2010,5))),colour=NA) +
    scale_fill_manual(values = palette5,
                      labels=substr(quintileBreaks(fishnetPopulation,"pop_2010"),1,4),
                      name="Quintile\nBreaks") +
    labs(title="Population, DVR: 2010",
         subtitle="Represented as fishnet gridcells; Boundaries omitted") +
    mapTheme(), ncol=2)

2.5 Highway distance

DVRHighways <-
  read_sf("DVR_highway2proj.shp") %>%
  st_transform(st_crs(DVR)) %>%
  st_intersection(DVR)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
ggplot() +
  geom_point(data=fishnet, 
             aes(x=xyC(fishnet)[,1], y=xyC(fishnet)[,2],colour=lc_change),size=1.5) +
  geom_sf(data=DVRHighways) +
  scale_colour_manual(values = palette2,na.translate = F,
                      labels=c("No Change","New Development")) +
  labs(title = "New Development and Highways",
       subtitle = "As fishnet centroids") +
  mapTheme()
## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x
## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

emptyRaster <- lc_change
emptyRaster[] <- NA


highway_raster <- 
  as(DVRHighways,'Spatial') %>%
  rasterize(.,emptyRaster)

plot(highway_raster)

highway_raster_distance <- distance(highway_raster)
names(highway_raster_distance) <- "distance_highways"

plot(highway_raster_distance)

highwayPoints <-
  rasterToPoints(highway_raster_distance) %>%
  as.data.frame() %>%
  st_as_sf(coords = c("x", "y"), crs = st_crs(DVR_fishnet))

plot(highwayPoints)

highwayPoints_fishnet <- 
  aggregate(highwayPoints, DVR_fishnet, mean) %>%
  mutate(distance_highways = ifelse(is.na(distance_highways),0,distance_highways))

ggplot() +
  geom_sf(data=DVR) +
  geom_point(data=highwayPoints_fishnet, aes(x=xyC(highwayPoints_fishnet)[,1], 
                                             y=xyC(highwayPoints_fishnet)[,2], 
                                             colour=factor(ntile(distance_highways,5))),size=1.5) +
  scale_colour_manual(values = palette5,
                      labels=substr(quintileBreaks(highwayPoints_fishnet,"distance_highways"),1,8),
                      name="Quintile\nBreaks") +
  geom_sf(data=DVRHighways, colour = "red") +
  labs(title = "Distance to Highways",
       subtitle = "As fishnet centroids; Highways visualized in red") +
  mapTheme()
## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

2.6 The spatial lag of development

nn_function <- function(measureFrom,measureTo,k) {
  #convert the sf layers to matrices
  measureFrom_Matrix <-
    as.matrix(measureFrom)
  measureTo_Matrix <-
    as.matrix(measureTo)
  nn <-   
    get.knnx(measureTo, measureFrom, k)$nn.dist
  output <-
    as.data.frame(nn) %>%
    rownames_to_column(var = "thisPoint") %>%
    gather(points, point_distance, V1:ncol(.)) %>%
    arrange(as.numeric(thisPoint)) %>%
    group_by(thisPoint) %>%
    summarize(pointDistance = mean(point_distance)) %>%
    arrange(as.numeric(thisPoint)) %>% 
    dplyr::select(-thisPoint) %>%
    pull()
  
  return(output)  
}



fishnet$lagDevelopment <-
  nn_function(xyC(fishnet),
              xyC(filter(aggregatedRasters,developed==1)),
              2)
## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x
ggplot() +
  geom_sf(data=DVR) +
  geom_point(data=fishnet, 
             aes(x=xyC(fishnet)[,1], y=xyC(fishnet)[,2], 
                 colour=factor(ntile(lagDevelopment,5))), size=1.5) +
  scale_colour_manual(values = palette5,
                      labels=substr(quintileBreaks(fishnet,"lagDevelopment"),1,7),
                      name="Quintile\nBreaks") +
  labs(title = "Spatial lag to 2001 development",
       subtitle = "As fishnet centroids") +
  mapTheme()
## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

# 2.7 MSA Counties

options(tigris_class = "sf")

PACounties <- 
  counties("Pennsylvania") %>%
  st_transform(st_crs(DVR)) %>%
  dplyr::select(NAME) %>%
  .[st_buffer(DVR,-1000), , op=st_intersects]
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |                                                                      |   1%
  |                                                                            
  |=                                                                     |   1%
  |                                                                            
  |=                                                                     |   2%
  |                                                                            
  |==                                                                    |   2%
  |                                                                            
  |==                                                                    |   3%
  |                                                                            
  |==                                                                    |   4%
  |                                                                            
  |===                                                                   |   4%
  |                                                                            
  |===                                                                   |   5%
  |                                                                            
  |====                                                                  |   5%
  |                                                                            
  |====                                                                  |   6%
  |                                                                            
  |=====                                                                 |   6%
  |                                                                            
  |=====                                                                 |   7%
  |                                                                            
  |=====                                                                 |   8%
  |                                                                            
  |======                                                                |   8%
  |                                                                            
  |======                                                                |   9%
  |                                                                            
  |=======                                                               |   9%
  |                                                                            
  |=======                                                               |  10%
  |                                                                            
  |=======                                                               |  11%
  |                                                                            
  |========                                                              |  11%
  |                                                                            
  |========                                                              |  12%
  |                                                                            
  |=========                                                             |  12%
  |                                                                            
  |=========                                                             |  13%
  |                                                                            
  |=========                                                             |  14%
  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |==========                                                            |  15%
  |                                                                            
  |===========                                                           |  15%
  |                                                                            
  |===========                                                           |  16%
  |                                                                            
  |============                                                          |  16%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |============                                                          |  18%
  |                                                                            
  |=============                                                         |  18%
  |                                                                            
  |=============                                                         |  19%
  |                                                                            
  |==============                                                        |  19%
  |                                                                            
  |==============                                                        |  20%
  |                                                                            
  |==============                                                        |  21%
  |                                                                            
  |===============                                                       |  21%
  |                                                                            
  |===============                                                       |  22%
  |                                                                            
  |================                                                      |  22%
  |                                                                            
  |================                                                      |  23%
  |                                                                            
  |================                                                      |  24%
  |                                                                            
  |=================                                                     |  24%
  |                                                                            
  |=================                                                     |  25%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |==================                                                    |  26%
  |                                                                            
  |===================                                                   |  26%
  |                                                                            
  |===================                                                   |  27%
  |                                                                            
  |===================                                                   |  28%
  |                                                                            
  |====================                                                  |  28%
  |                                                                            
  |====================                                                  |  29%
  |                                                                            
  |=====================                                                 |  29%
  |                                                                            
  |=====================                                                 |  30%
  |                                                                            
  |=====================                                                 |  31%
  |                                                                            
  |======================                                                |  31%
  |                                                                            
  |======================                                                |  32%
  |                                                                            
  |=======================                                               |  32%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |=======================                                               |  34%
  |                                                                            
  |========================                                              |  34%
  |                                                                            
  |========================                                              |  35%
  |                                                                            
  |=========================                                             |  35%
  |                                                                            
  |=========================                                             |  36%
  |                                                                            
  |==========================                                            |  36%
  |                                                                            
  |==========================                                            |  37%
  |                                                                            
  |==========================                                            |  38%
  |                                                                            
  |===========================                                           |  38%
  |                                                                            
  |===========================                                           |  39%
  |                                                                            
  |============================                                          |  39%
  |                                                                            
  |============================                                          |  40%
  |                                                                            
  |============================                                          |  41%
  |                                                                            
  |=============================                                         |  41%
  |                                                                            
  |=============================                                         |  42%
  |                                                                            
  |==============================                                        |  42%
  |                                                                            
  |==============================                                        |  43%
  |                                                                            
  |==============================                                        |  44%
  |                                                                            
  |===============================                                       |  44%
  |                                                                            
  |===============================                                       |  45%
  |                                                                            
  |================================                                      |  45%
  |                                                                            
  |================================                                      |  46%
  |                                                                            
  |=================================                                     |  46%
  |                                                                            
  |=================================                                     |  47%
  |                                                                            
  |=================================                                     |  48%
  |                                                                            
  |==================================                                    |  48%
  |                                                                            
  |==================================                                    |  49%
  |                                                                            
  |===================================                                   |  49%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |===================================                                   |  51%
  |                                                                            
  |====================================                                  |  51%
  |                                                                            
  |====================================                                  |  52%
  |                                                                            
  |=====================================                                 |  52%
  |                                                                            
  |=====================================                                 |  53%
  |                                                                            
  |=====================================                                 |  54%
  |                                                                            
  |======================================                                |  54%
  |                                                                            
  |======================================                                |  55%
  |                                                                            
  |=======================================                               |  55%
  |                                                                            
  |=======================================                               |  56%
  |                                                                            
  |========================================                              |  56%
  |                                                                            
  |========================================                              |  57%
  |                                                                            
  |========================================                              |  58%
  |                                                                            
  |=========================================                             |  58%
  |                                                                            
  |=========================================                             |  59%
  |                                                                            
  |==========================================                            |  59%
  |                                                                            
  |==========================================                            |  60%
  |                                                                            
  |==========================================                            |  61%
  |                                                                            
  |===========================================                           |  61%
  |                                                                            
  |===========================================                           |  62%
  |                                                                            
  |============================================                          |  62%
  |                                                                            
  |============================================                          |  63%
  |                                                                            
  |============================================                          |  64%
  |                                                                            
  |=============================================                         |  64%
  |                                                                            
  |=============================================                         |  65%
  |                                                                            
  |==============================================                        |  65%
  |                                                                            
  |==============================================                        |  66%
  |                                                                            
  |===============================================                       |  66%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |===============================================                       |  68%
  |                                                                            
  |================================================                      |  68%
  |                                                                            
  |================================================                      |  69%
  |                                                                            
  |=================================================                     |  69%
  |                                                                            
  |=================================================                     |  70%
  |                                                                            
  |=================================================                     |  71%
  |                                                                            
  |==================================================                    |  71%
  |                                                                            
  |==================================================                    |  72%
  |                                                                            
  |===================================================                   |  72%
  |                                                                            
  |===================================================                   |  73%
  |                                                                            
  |===================================================                   |  74%
  |                                                                            
  |====================================================                  |  74%
  |                                                                            
  |====================================================                  |  75%
  |                                                                            
  |=====================================================                 |  75%
  |                                                                            
  |=====================================================                 |  76%
  |                                                                            
  |======================================================                |  76%
  |                                                                            
  |======================================================                |  77%
  |                                                                            
  |======================================================                |  78%
  |                                                                            
  |=======================================================               |  78%
  |                                                                            
  |=======================================================               |  79%
  |                                                                            
  |========================================================              |  79%
  |                                                                            
  |========================================================              |  80%
  |                                                                            
  |========================================================              |  81%
  |                                                                            
  |=========================================================             |  81%
  |                                                                            
  |=========================================================             |  82%
  |                                                                            
  |==========================================================            |  82%
  |                                                                            
  |==========================================================            |  83%
  |                                                                            
  |==========================================================            |  84%
  |                                                                            
  |===========================================================           |  84%
  |                                                                            
  |===========================================================           |  85%
  |                                                                            
  |============================================================          |  85%
  |                                                                            
  |============================================================          |  86%
  |                                                                            
  |=============================================================         |  86%
  |                                                                            
  |=============================================================         |  87%
  |                                                                            
  |=============================================================         |  88%
  |                                                                            
  |==============================================================        |  88%
  |                                                                            
  |==============================================================        |  89%
  |                                                                            
  |===============================================================       |  89%
  |                                                                            
  |===============================================================       |  90%
  |                                                                            
  |===============================================================       |  91%
  |                                                                            
  |================================================================      |  91%
  |                                                                            
  |================================================================      |  92%
  |                                                                            
  |=================================================================     |  92%
  |                                                                            
  |=================================================================     |  93%
  |                                                                            
  |=================================================================     |  94%
  |                                                                            
  |==================================================================    |  94%
  |                                                                            
  |==================================================================    |  95%
  |                                                                            
  |===================================================================   |  95%
  |                                                                            
  |===================================================================   |  96%
  |                                                                            
  |====================================================================  |  96%
  |                                                                            
  |====================================================================  |  97%
  |                                                                            
  |====================================================================  |  98%
  |                                                                            
  |===================================================================== |  98%
  |                                                                            
  |===================================================================== |  99%
  |                                                                            
  |======================================================================|  99%
  |                                                                            
  |======================================================================| 100%
NjCounties <- 
  counties("New Jersey") %>%
  st_transform(st_crs(DVR)) %>%
  dplyr::select(NAME) %>%
  .[st_buffer(DVR,-1000), , op=st_intersects]

studyAreaCounties <- rbind(PACounties, NjCounties)
  


ggplot() +
  geom_sf(data=studyAreaCounties) +
  labs(title = "Study area counties") +
  mapTheme()

3. Exploratory Analysis

head(fishnet)
## Simple feature collection with 6 features and 2 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 2738314 ymin: 78345.53 xmax: 2754314 ymax: 86345.53
## CRS:            2272
##   lc_change                       geometry lagDevelopment
## 1         0 POLYGON ((2742314 78345.53,...      10000.000
## 2         0 POLYGON ((2746314 78345.53,...      10796.691
## 3         0 POLYGON ((2750314 78345.53,...      12867.957
## 4         0 POLYGON ((2738314 82345.53,...       7300.563
## 5         0 POLYGON ((2742314 82345.53,...       6000.000
## 6         0 POLYGON ((2746314 82345.53,...       7300.563
dat <- 
  cbind(
    fishnet, highwayPoints_fishnet, fishnetPopulation, aggregatedRasters) %>%
  dplyr::select(lc_change, developed, forest, farm, wetlands, otherUndeveloped, water,
                pop_2000, pop_2010, pop_Change, distance_highways,lagDevelopment) %>%
  st_join(studyAreaCounties) %>%
  mutate(developed10 = ifelse(lc_change == 1 & developed == 1, 0, developed)) %>%
  filter(water == 0) 

dat %>%
  mutate(distance_highways = replace_na(distance_highways, 0)) 
## Simple feature collection with 7182 features and 14 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 2422314 ymin: 78345.53 xmax: 2918314 ymax: 478345.5
## CRS:            2272
## First 10 features:
##    lc_change developed forest farm wetlands otherUndeveloped water  pop_2000
## 1          0         0      0    0        0                0     0  31.47190
## 2          0         0      1    0        0                0     0 138.88617
## 3          0         0      0    1        0                0     0  15.25962
## 4          0         0      0    0        0                0     0  51.22063
## 5          0         0      1    0        0                0     0 174.92009
## 6          0         0      0    1        0                0     0 185.58596
## 7          0         0      1    0        0                0     0 141.19596
## 8          0         0      0    0        0                0     0   4.74980
## 9          0         0      0    0        0                0     0 217.97177
## 10         0         0      0    0        0                0     0 525.34027
##      pop_2010 pop_Change distance_highways lagDevelopment       NAME
## 1   35.679977   4.208077          26731.46      10000.000 Gloucester
## 2  159.361500  20.475328          30247.79      10796.691 Gloucester
## 3   18.054369   2.794750          32808.33      12867.957 Gloucester
## 4   58.230532   7.009906          19783.39       7300.563 Gloucester
## 5  199.249651  24.329560          23372.28       6000.000 Gloucester
## 6  211.646562  26.060602          27336.86       7300.563 Gloucester
## 7  161.818594  20.622635          30147.99      10128.990 Gloucester
## 8    5.313537   0.563737          33045.35      13535.658 Gloucester
## 9  213.771697  -4.200072          13919.40       6828.427 Gloucester
## 10 502.601443 -22.738826          16404.17       4828.427 Gloucester
##                          geometry developed10
## 1  POLYGON ((2742314 78345.53,...           1
## 2  POLYGON ((2746314 78345.53,...           1
## 3  POLYGON ((2750314 78345.53,...           1
## 4  POLYGON ((2738314 82345.53,...           1
## 5  POLYGON ((2742314 82345.53,...           1
## 6  POLYGON ((2746314 82345.53,...           1
## 7  POLYGON ((2750314 82345.53,...           1
## 8  POLYGON ((2754314 82345.53,...           1
## 9  POLYGON ((2734314 86345.53,...           1
## 10 POLYGON ((2738314 86345.53,...           1
dat %>%
  mutate(lc_change = replace_na(lc_change, 0)) 
## Simple feature collection with 7182 features and 14 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 2422314 ymin: 78345.53 xmax: 2918314 ymax: 478345.5
## CRS:            2272
## First 10 features:
##    lc_change developed forest farm wetlands otherUndeveloped water  pop_2000
## 1          0         0      0    0        0                0     0  31.47190
## 2          0         0      1    0        0                0     0 138.88617
## 3          0         0      0    1        0                0     0  15.25962
## 4          0         0      0    0        0                0     0  51.22063
## 5          0         0      1    0        0                0     0 174.92009
## 6          0         0      0    1        0                0     0 185.58596
## 7          0         0      1    0        0                0     0 141.19596
## 8          0         0      0    0        0                0     0   4.74980
## 9          0         0      0    0        0                0     0 217.97177
## 10         0         0      0    0        0                0     0 525.34027
##      pop_2010 pop_Change distance_highways lagDevelopment       NAME
## 1   35.679977   4.208077          26731.46      10000.000 Gloucester
## 2  159.361500  20.475328          30247.79      10796.691 Gloucester
## 3   18.054369   2.794750          32808.33      12867.957 Gloucester
## 4   58.230532   7.009906          19783.39       7300.563 Gloucester
## 5  199.249651  24.329560          23372.28       6000.000 Gloucester
## 6  211.646562  26.060602          27336.86       7300.563 Gloucester
## 7  161.818594  20.622635          30147.99      10128.990 Gloucester
## 8    5.313537   0.563737          33045.35      13535.658 Gloucester
## 9  213.771697  -4.200072          13919.40       6828.427 Gloucester
## 10 502.601443 -22.738826          16404.17       4828.427 Gloucester
##                          geometry developed10
## 1  POLYGON ((2742314 78345.53,...           1
## 2  POLYGON ((2746314 78345.53,...           1
## 3  POLYGON ((2750314 78345.53,...           1
## 4  POLYGON ((2738314 82345.53,...           1
## 5  POLYGON ((2742314 82345.53,...           1
## 6  POLYGON ((2746314 82345.53,...           1
## 7  POLYGON ((2750314 82345.53,...           1
## 8  POLYGON ((2754314 82345.53,...           1
## 9  POLYGON ((2734314 86345.53,...           1
## 10 POLYGON ((2738314 86345.53,...           1
dat %>%
  dplyr::select(distance_highways,lagDevelopment,lc_change) %>%
  gather(Variable, Value, -lc_change, -geometry) %>%
    ggplot(., aes(lc_change, Value, fill=lc_change,), Variable[!is.na(Variable$lc_change),]) + 
  geom_bar(position = "dodge", stat = "summary", fun.y = "mean") +
  facet_wrap(~Variable) +
  scale_fill_manual(values = palette2,na.translate = F,
                    labels=c("No Change","New Development"),
                    name="") +
  labs(title="New development as a function of the countinuous variables") +
  plotTheme() 
## Warning: Ignoring unknown parameters: fun.y
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`

dat %>%
  dplyr::select(pop_2000,pop_2010,pop_Change,lc_change) %>%
  gather(Variable, Value, -lc_change, -geometry) %>%
  ggplot(., aes(lc_change, Value, fill=lc_change)) + 
  geom_bar(position = "dodge", stat = "summary", fun.y = "mean") +
  facet_wrap(~Variable) +
  scale_fill_manual(values = palette2,na.translate = F,
                    labels=c("No Change","New Development"),
                    name="") +
  labs(title="New development as a function of factor variables") +
  plotTheme()
## Warning: Ignoring unknown parameters: fun.y
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`

dat %>%
  dplyr::select(lc_change:otherUndeveloped,developed) %>%
  gather(Land_Cover_Type, Value, -lc_change, -geometry) %>%
  st_set_geometry(NULL) %>%
  group_by(lc_change, Land_Cover_Type) %>%
  summarize(n = sum(as.numeric(Value))) %>%
  ungroup() %>%
  mutate(Conversion_Rate = paste0(round(100 * n/sum(n), 2), "%")) %>%
  filter(lc_change == 1) %>%
  dplyr::select(Land_Cover_Type,Conversion_Rate) %>%
  kable() %>% kable_styling(full_width = F)
Land_Cover_Type Conversion_Rate
developed 1.07%
farm 1.16%
forest 0.67%
otherUndeveloped 0.05%
wetlands 0.26%
  1. Predicting for 2010
set.seed(3456)
trainIndex <- 
  createDataPartition(dat$developed, p = .50,
                      list = FALSE,
                      times = 1)
datTrain <- dat[ trainIndex,]
datTest  <- dat[-trainIndex,]

nrow(dat)
## [1] 7182
Model1 <- glm(lc_change ~ wetlands + forest  + farm + otherUndeveloped, 
              family="binomial"(link="logit"), data = datTrain)

Model2 <- glm(lc_change ~ wetlands + forest  + farm + otherUndeveloped + lagDevelopment, 
              family="binomial"(link="logit"), data = datTrain)

Model3 <- glm(lc_change ~ wetlands + forest  + farm + otherUndeveloped + lagDevelopment + pop_2000, 
              family="binomial"(link="logit"), data = datTrain)          

Model4 <- glm(lc_change ~ wetlands + forest  + farm + otherUndeveloped + lagDevelopment + pop_2000 + 
                pop_2010, 
              family="binomial"(link="logit"), data = datTrain) 


Model5 <- glm(lc_change ~ wetlands + forest  + farm + otherUndeveloped + lagDevelopment + pop_Change, 
              family="binomial"(link="logit"), data = datTrain)              

Model6 <- glm(lc_change ~ wetlands + forest  + farm + otherUndeveloped + lagDevelopment + pop_Change + 
                distance_highways, 
              family="binomial"(link="logit"), data = datTrain) 


modelList <- paste0("Model", 1:6)
map_dfc(modelList, function(x)pR2(get(x)))[4,] %>%
  setNames(paste0("Model",1:6)) %>%
  gather(Model,McFadden) %>%
  ggplot(aes(Model,McFadden)) +
  geom_bar(stat="identity") +
  labs(title= "McFadden R-Squared by Model") +
  plotTheme()
## fitting null model for pseudo-r2
## fitting null model for pseudo-r2
## fitting null model for pseudo-r2
## fitting null model for pseudo-r2
## fitting null model for pseudo-r2
## fitting null model for pseudo-r2

testSetProbs <- 
  data.frame(class = datTest$lc_change,
             probs = predict(Model6, datTest, type="response")) 

ggplot(testSetProbs, aes(probs)) +
  geom_density(aes(fill=class), alpha=0.5) +
  scale_fill_manual(values = palette2,na.translate = F,
                    labels=c("No Change","New Development")) +
  labs(title = "Histogram of test set predicted probabilities",
       x="Predicted Probabilities",y="Density") +
  plotTheme()

head(testSetProbs)
##    class      probs
## 3      0 0.02294882
## 5      0 0.01252229
## 6      0 0.04237250
## 9      0 0.01594626
## 10     0 0.01863234
## 11     0 0.02471212
options(yardstick.event_first = FALSE)

testSetProbs <- 
  testSetProbs %>% 
  mutate(predClass_05 = as.factor(ifelse(testSetProbs$probs >= 0.05 ,1,0)),
         predClass_17 = as.factor(ifelse(testSetProbs$probs >= 0.17 ,1,0))) 

testSetProbs %>%
  dplyr::select(-probs) %>%
  gather(Variable, Value, -class) %>%
  group_by(Variable) %>%
  summarize(Sensitivity = round(yardstick::sens_vec(class,factor(Value)),2),
            Specificity = round(yardstick::spec_vec(class,factor(Value)),2),
            Accuracy = round(yardstick::accuracy_vec(class,factor(Value)),2)) %>% 
  kable() %>%
  kable_styling(full_width = F)
Variable Sensitivity Specificity Accuracy
predClass_05 0.37 0.88 0.86
predClass_17 0.00 1.00 0.97
predsForMap <-         
  dat %>%
  mutate(probs = predict(Model6, dat, type="response") ,
         Threshold_5_Pct = as.factor(ifelse(probs >= 0.05 ,1,0)),
         Threshold_17_Pct =  as.factor(ifelse(probs >= 0.17 ,1,0))) %>%
  dplyr::select(lc_change,Threshold_5_Pct,Threshold_17_Pct) %>%
  gather(Variable,Value, -geometry) %>%
  st_cast("POLYGON")
## Warning: attributes are not identical across measure variables;
## they will be dropped
ggplot() +
  geom_point(data=predsForMap, aes(x=xyC(predsForMap)[,1], y=xyC(predsForMap)[,2], colour=Value)) +
  facet_wrap(~Variable) +
  scale_colour_manual(values = palette2, labels=c("No Change","New Development"),na.translate = F,
                      name="") +
  labs(title="Development predictions - Low threshold") + mapTheme()
## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x
## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

ConfusionMatrix.metrics <-
  dat %>%
  mutate(probs = predict(Model6, dat, type="response") ,
         Threshold_5_Pct = as.factor(ifelse(probs >= 0.05 ,1,0)),
         Threshold_17_Pct =  as.factor(ifelse(probs >= 0.17 ,1,0))) %>%
  mutate(TrueP_05 = ifelse(lc_change  == 1 & Threshold_5_Pct == 1, 1,0),
         TrueN_05 = ifelse(lc_change  == 0 & Threshold_5_Pct == 0, 1,0),
         TrueP_17 = ifelse(lc_change  == 1 & Threshold_17_Pct == 1, 1,0),
         TrueN_17 = ifelse(lc_change  == 0 & Threshold_17_Pct == 0, 1,0)) %>%
  dplyr::select(., starts_with("True")) %>%
  gather(Variable, Value, -geometry) %>%
  st_cast("POLYGON") 

ggplot(data=ConfusionMatrix.metrics) +
  geom_point(aes(x=xyC(ConfusionMatrix.metrics)[,1], 
                 y=xyC(ConfusionMatrix.metrics)[,2], colour = as.factor(Value))) +
  facet_wrap(~Variable) +
  scale_colour_manual(values = palette2, labels=c("Correct","Incorrect"),na.translate = F,
                      name="") +
  labs(title="Development predictions - Low threshold") + mapTheme()
## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

4.4 Generalizability

spatialCV <- function(dataFrame, uniqueID, dependentVariable, modelName) {
  
  #initialize a data frame 
  endList <- list()
  
  #create a list that is all the spatial group unqiue ids in the data frame (ie counties)    
  uniqueID_List <- unique(dataFrame[[uniqueID]])  
  x <- 1
  y <- length(uniqueID_List)
  
  #create a counter and while it is less than the number of counties...  
  while(x <= y) 
  {
    #call a current county    
    currentUniqueID <- uniqueID_List[x]
    #create a training set comprised of units not in that county and a test set of units
    #that are that county
    training <- dataFrame[ which(dataFrame[[uniqueID]] != uniqueID_List[x]),]
    testing <- dataFrame[ which(dataFrame[[uniqueID]] == uniqueID_List[x]),]
    #create seperate xy vectors
    trainingX <- training[ , -which(names(training) %in% c(dependentVariable))]
    testingX <- testing[ , -which(names(testing) %in% c(dependentVariable))]
    
    trainY <- training[[dependentVariable]]
    testY <- testing[[dependentVariable]]
    #Calculate predictions on the test county as part of a data frame including the observed
    #outcome and the unique county ID    
    thisPrediction <- 
      data.frame(class = testY,
                 probs = predict(modelName, testingX, type="response"),
                 county = currentUniqueID) 
    
    #Row bind the predictions to a data farme
    endList <- rbind(endList, thisPrediction)
    #iterate counter    
    x <- x + 1 
  } 
  #return the final list of counties and associated predictions  
  return (as.data.frame(endList))
}



spatialCV_counties <-
  spatialCV(dat,"NAME","lc_change", Model6) %>%
  mutate(predClass = as.factor(ifelse(probs >= 0.17 ,1,0)))


spatialCV_metrics <-
  spatialCV_counties %>% 
  group_by(county) %>% 
  summarize(Observed_Change = sum(as.numeric(as.character(class))),
            Sensitivity = round(yardstick::sens_vec(class,predClass),2),
            Specificity = round(yardstick::spec_vec(class,predClass),2),
            Accuracy = round(yardstick::accuracy_vec(class,predClass),2)) 

spatialCV_metrics %>%
  kable() %>%
  kable_styling(full_width = F)
county Observed_Change Sensitivity Specificity Accuracy
Bucks 24 0 1.00 0.98
Burlington 47 0 1.00 0.97
Camden 7 0 1.00 0.98
Chester 35 0 1.00 0.98
Delaware 9 0 1.00 0.98
Gloucester 36 0 1.00 0.94
Mercer 17 0 1.00 0.96
Montgomery 37 0 1.00 0.96
Philadelphia 5 0 0.98 0.96
  1. Prediction of land cover change in 2020 Demand-side prediction of development
dat <-
  dat %>%
  mutate(lagDevelopment = nn_function(xyC(.), xyC(filter(.,developed10 == 1)),2))
## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x
countyPopulation_2020 <- 
  data.frame(
    NAME = 
      c("Bucks","Montgomery","Chester","Philadelphia","Delaware",
        "Mercer","Camden","Gloucester","Burlington"),
    county_projection_2020 = 
      c(640495, 840934, 543702, 1594787, 568337, 377328 , 514006, 307766, 459344)) %>%
  left_join(
    dat %>%
      st_set_geometry(NULL) %>%
      group_by(NAME) %>%
      summarize(county_population_2010 = round(sum(pop_2010))))
## Joining, by = "NAME"
countyPopulation_2020 %>%
  gather(Variable,Value, -NAME) %>%
  ggplot(aes(reorder(NAME,-Value),Value)) +
  geom_bar(aes(fill=Variable), stat = "identity") +
  scale_fill_manual(values = palette2,
                    labels=c("2020","2010"),
                    name="Population") +
  labs(title="Population Change by county: 2010 - 2020",
       x="County", y="Population") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  plotTheme()

dat_infill <-
  dat %>%
  #calculate population change
  left_join(countyPopulation_2020) %>%
  mutate(proportion_of_county_pop = pop_2010 / county_population_2010,
         pop_2020.infill = proportion_of_county_pop * county_projection_2020,
         pop_Change = round(pop_2020.infill - pop_2010),2) %>%
  dplyr::select(-county_projection_2020, -county_population_2010, 
                -proportion_of_county_pop, -pop_2020.infill) %>%
  #predict for 2020
  mutate(predict_2020.infill = predict(Model6,. , type="response"))
## Joining, by = "NAME"
dat_infill %>%
  ggplot() +  
  geom_point(aes(x=xyC(dat_infill)[,1], y=xyC(dat_infill)[,2], colour = factor(ntile(predict_2020.infill,5)))) +
  scale_colour_manual(values = palette5,
                      labels=substr(quintileBreaks(dat_infill,"predict_2020.infill"),1,4),
                      name="Quintile\nBreaks") +
  geom_sf(data=studyAreaCounties, fill=NA, colour="black", size=1.5) +
  labs(title= "Development Demand in 2020: Predicted Probabilities") +
  mapTheme()
## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

6. Comparing predicted development demand and environmental sensitivity

lc_2011 <- raster("lc_2011_1kproj.tif")

developed11 <- lc_2011 == 21 | lc_2011 == 22 | lc_2011 == 23 | lc_2011 == 24
forest11 <- lc_2011 == 41 | lc_2011 == 42 | lc_2011 == 43 
farm11 <- lc_2011 == 81 | lc_2011 == 82 
wetlands11 <- lc_2011 == 90 | lc_2011 == 95 
otherUndeveloped11 <- lc_2011 == 52 | lc_2011 == 71 | lc_2011 == 31 
water11 <- lc_2011 == 11

names(developed11) <- "developed11"
names(forest11) <- "forest11"
names(farm11) <- "farm11"
names(wetlands11) <- "wetlands11"
names(otherUndeveloped11) <- "otherUndeveloped11"
names(water11) <- "water11"

ggplot() +
  geom_sf(data=DVR) +
  geom_raster(data = rbind(rast(lc_2001) %>% mutate(label = "2001"),
                           rast(lc_2011) %>% mutate(label = "2011")) %>% 
                na.omit %>% filter(value > 0), 
              aes(x,y,fill=as.factor(value))) +
  facet_wrap(~label) +
  scale_fill_viridis(discrete=TRUE, name ="") +
  labs(title = "Land Cover, 2001 & 2011") +
  mapTheme() + theme(legend.position = "none")

theRasterList11 <- c(developed11,forest11,farm11,wetlands11,otherUndeveloped11,water11)

dat2 <-
  aggregateRaster(theRasterList11, dat) %>%
  dplyr::select(developed11,forest11,farm11,wetlands11,otherUndeveloped11,water11) %>%
  st_set_geometry(NULL) %>%
  bind_cols(.,dat) %>%
  st_sf() %>%
  st_cast("POLYGON")

dat2 %>%
  gather(var,value,developed11:water11) %>%
  st_centroid() %>%
  mutate(X = st_coordinates(.)[,1],
         Y = st_coordinates(.)[,2]) %>%
  ggplot() +
  geom_sf(data=DVR) +
  geom_point(aes(X,Y, colour=as.factor(value))) +
  facet_wrap(~var) +
  scale_colour_manual(values = palette2,
                      labels=c("Other","Land Cover"),
                      name = "") +
  labs(title = "Land cover types, 2011",
       subtitle = "As fishnet centroids") +
  mapTheme()
## Warning in st_centroid.sf(.): st_centroid assumes attributes are constant over
## geometries of x

# 6.3

dat2 <-
  dat2 %>%
  mutate(sensitive_lost11 = ifelse(forest == 1 & forest11 == 0 |
                                     wetlands == 1 & wetlands11 == 0,1,0))

ggplot() +
  geom_point(data=dat2, aes(x=xyC(dat2)[,1], y=xyC(dat2)[,2], colour=as.factor(sensitive_lost11))) +
  scale_colour_manual(values = palette2,
                      labels=c("No Change","Sensitive Lost"),
                      name = "") +
  labs(title = "Sensitive lands lost: 2001 - 2011",
       subtitle = "As fishnet centroids") +
  mapTheme()
## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

# 6.4

#install.packages("igraph")
#library(igraph)
sensitiveRegions <- 
  clump(wetlands11 + forest11) %>%
  rasterToPolygons() %>%
  st_as_sf() %>%
  group_by(clumps) %>% summarize() %>%
  mutate(Acres = as.numeric(st_area(.) * 0.0000229568)) %>%
  filter(Acres > 3954)  %>%
  dplyr::select() %>%
  rasterize(.,emptyRaster) 
## Loading required namespace: igraph
sensitiveRegions[sensitiveRegions > 0] <- 1  
names(sensitiveRegions) <- "sensitiveRegions"
dat2 <-
  aggregateRaster(c(sensitiveRegions), dat2) %>%
  dplyr::select(sensitiveRegions) %>%
  st_set_geometry(NULL) %>%
  bind_cols(.,dat2) %>%
  st_sf()

ggplot() +
  geom_point(data=dat2, aes(x=xyC(dat2)[,1], y=xyC(dat2)[,2], colour=as.factor(sensitiveRegions))) +
  scale_colour_manual(values = palette2,
                      labels=c("Other","Sensitive Regions"),
                      name="") +
  labs(title = "Sensitive regions",
       subtitle = "Continous areas of either wetlands or forests\ngreater than 1 acre") +
  mapTheme()
## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

county_specific_metrics <- 
  dat2 %>%
  #predict development demand from our model
  mutate(Development_Demand = predict(Model6, dat2, type="response")) %>%
  #get a count count of grid cells by county which we can use to calculate rates below
  left_join(st_set_geometry(dat, NULL) %>% group_by(NAME) %>% summarize(count = n())) %>%
  #calculate summary statistics by county
  group_by(NAME) %>%
  summarize(Total_Farmland = sum(farm11) / max(count),
            Total_Forest = sum(forest11) / max(count),
            Total_Wetlands = sum(wetlands11) / max(count),
            Total_Undeveloped = sum(otherUndeveloped11) / max(count),
            Sensitive_Land_Lost = sum(sensitive_lost11) / max(count),
            Sensitive_Regions = sum(sensitiveRegions) / max(count),
            Mean_Development_Demand = mean(Development_Demand)) %>%
  #get population data by county
  left_join(countyPopulation_2020 %>% 
              mutate(Population_Change = county_projection_2020 - county_population_2010,
                     Population_Change_Rate = Population_Change / county_projection_2020) %>%
              dplyr::select(NAME,Population_Change_Rate))
## Joining, by = "NAME"
## Joining, by = "NAME"
county_specific_metrics %>%
  gather(Variable, Value, -NAME, -geometry) %>%
  mutate(Variable = factor(Variable, levels=c("Population_Change_Rate","Mean_Development_Demand",
                                              "Total_Farmland","Total_Undeveloped","Total_Forest",
                                              "Total_Wetlands","Sensitive_Land_Lost","Sensitive_Regions",
                                              ordered = TRUE))) %>%
  mutate(Planning_Designation = case_when(
    Variable == "Population_Change_Rate" | Variable == "Mean_Development_Demand" ~ "Demand-Side",
    Variable == "Total_Farmland" | Variable == "Total_Undeveloped"               ~ "Suitable",
    TRUE                                                                         ~ "Not Suitable")) %>%
  ggplot(aes(x=Variable, y=Value, fill=Planning_Designation)) +
  geom_bar(stat="identity", position=position_dodge(), colour="black") +
  facet_wrap(~NAME, ncol=5) +
  coord_flip() +
  scale_y_continuous(breaks = seq(.25, 1, by = .25)) +
  geom_vline(xintercept = 2.5) + geom_vline(xintercept = 4.5) +
  scale_fill_manual(values=c("black","red","darkgreen")) +
  labs(title= "County Specific Allocation Metrics", subtitle= "As rates", x="Indicator", y="Rate") +
  plotTheme() + theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position="bottom")

7. Allocation

# Philadelphia 2020
Philadelphia <-
  dat2 %>%
  mutate(Development_Demand = predict(Model6, dat2, type="response")) %>%
  filter(NAME == "Philadelphia") 

Philadelphia_landUse <- rbind(
  filter(Philadelphia, forest11 == 1 | wetlands11 == 1 ) %>%
    dplyr::select() %>% mutate(Land_Use = "Not Suitable"),
  filter(Philadelphia, developed11 == 1) %>%
    dplyr::select() %>% mutate(Land_Use = "Developed"))

grid.arrange(
  ggplot() +
    geom_sf(data=Philadelphia, aes(fill=factor(ntile(Development_Demand,5))), colour=NA) +
    geom_point(data=Philadelphia_landUse, aes(x=xyC(Philadelphia_landUse)[,1], 
                                          y=xyC(Philadelphia_landUse)[,2], colour=Land_Use),
               shape = 15, size = 2) +
    geom_sf(data=st_intersection(DVRHighways,filter(studyAreaCounties, NAME=="Philadelphia")), size=2) +
    scale_fill_manual(values = palette5, name="Development\nDemand",
                      labels=substr(quintileBreaks(Philadelphia,"Development_Demand"),1,5)) +
    scale_colour_manual(values = c("black","red")) + 
    labs(title = "Development Potential, 2020: Philadelphia") + mapTheme() +
    guides(fill = guide_legend(order = 1), colour = guide_legend(order = 2)),
  
  ggplot() +
    geom_sf(data=Philadelphia, aes(fill=factor(ntile(pop_Change,5))), colour=NA) +
    geom_point(data=Philadelphia_landUse, aes(x=xyC(Philadelphia_landUse)[,1], 
                                          y=xyC(Philadelphia_landUse)[,2], colour=Land_Use),
               shape = 15, size = 2) +
    geom_sf(data=st_intersection(DVRHighways,filter(studyAreaCounties, NAME=="Philadelphia")), size=2) +
    scale_fill_manual(values = palette5, name="Population\nChange",
                      labels=substr(quintileBreaks(Philadelphia,"pop_Change"),1,5)) +
    scale_colour_manual(values = c("black","red")) + 
    labs(title = "Projected Population, 2020: Philadelphia") + mapTheme() +
    guides(fill = guide_legend(order = 1), colour = guide_legend(order = 2)), ncol=2)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

# Bucks
Philadelphia <-
  dat2 %>%
  mutate(Development_Demand = predict(Model6, dat2, type="response")) %>%
  filter(NAME == "Bucks") 

Philadelphia_landUse <- rbind(
  filter(Philadelphia, forest11 == 1 | wetlands11 == 1 ) %>%
    dplyr::select() %>% mutate(Land_Use = "Not Suitable"),
  filter(Philadelphia, developed11 == 1) %>%
    dplyr::select() %>% mutate(Land_Use = "Developed"))

grid.arrange(
  ggplot() +
    geom_sf(data=Philadelphia, aes(fill=factor(ntile(Development_Demand,5))), colour=NA) +
    geom_point(data=Philadelphia_landUse, aes(x=xyC(Philadelphia_landUse)[,1], 
                                              y=xyC(Philadelphia_landUse)[,2], colour=Land_Use),
               shape = 15, size = 2) +
    geom_sf(data=st_intersection(DVRHighways,filter(studyAreaCounties, NAME=="Bucks")), size=2) +
    scale_fill_manual(values = palette5, name="Development\nDemand",
                      labels=substr(quintileBreaks(Philadelphia,"Development_Demand"),1,5)) +
    scale_colour_manual(values = c("black","red")) + 
    labs(title = "Development Potential, 2020: Bucks") + mapTheme() +
    guides(fill = guide_legend(order = 1), colour = guide_legend(order = 2)),
  
  ggplot() +
    geom_sf(data=Philadelphia, aes(fill=factor(ntile(pop_Change,5))), colour=NA) +
    geom_point(data=Philadelphia_landUse, aes(x=xyC(Philadelphia_landUse)[,1], 
                                              y=xyC(Philadelphia_landUse)[,2], colour=Land_Use),
               shape = 15, size = 2) +
    geom_sf(data=st_intersection(DVRHighways,filter(studyAreaCounties, NAME=="Bucks")), size=2) +
    scale_fill_manual(values = palette5, name="Population\nChange",
                      labels=substr(quintileBreaks(Philadelphia,"pop_Change"),1,5)) +
    scale_colour_manual(values = c("black","red")) + 
    labs(title = "Projected Population, 2020: Bucks") + mapTheme() +
    guides(fill = guide_legend(order = 1), colour = guide_legend(order = 2)), ncol=2)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

# Montgomery 2020
Philadelphia <-
  dat2 %>%
  mutate(Development_Demand = predict(Model6, dat2, type="response")) %>%
  filter(NAME == "Montgomery") 

Philadelphia_landUse <- rbind(
  filter(Philadelphia, forest11 == 1 | wetlands11 == 1 ) %>%
    dplyr::select() %>% mutate(Land_Use = "Not Suitable"),
  filter(Philadelphia, developed11 == 1) %>%
    dplyr::select() %>% mutate(Land_Use = "Developed"))

grid.arrange(
  ggplot() +
    geom_sf(data=Philadelphia, aes(fill=factor(ntile(Development_Demand,5))), colour=NA) +
    geom_point(data=Philadelphia_landUse, aes(x=xyC(Philadelphia_landUse)[,1], 
                                              y=xyC(Philadelphia_landUse)[,2], colour=Land_Use),
               shape = 15, size = 2) +
    geom_sf(data=st_intersection(DVRHighways,filter(studyAreaCounties, NAME=="Montgomery")), size=2) +
    scale_fill_manual(values = palette5, name="Development\nDemand",
                      labels=substr(quintileBreaks(Philadelphia,"Development_Demand"),1,5)) +
    scale_colour_manual(values = c("black","red")) + 
    labs(title = "Development Potential, 2020: Montgomery") + mapTheme() +
    guides(fill = guide_legend(order = 1), colour = guide_legend(order = 2)),
  
  ggplot() +
    geom_sf(data=Philadelphia, aes(fill=factor(ntile(pop_Change,5))), colour=NA) +
    geom_point(data=Philadelphia_landUse, aes(x=xyC(Philadelphia_landUse)[,1], 
                                              y=xyC(Philadelphia_landUse)[,2], colour=Land_Use),
               shape = 15, size = 2) +
    geom_sf(data=st_intersection(DVRHighways,filter(studyAreaCounties, NAME=="Montgomery")), size=2) +
    scale_fill_manual(values = palette5, name="Population\nChange",
                      labels=substr(quintileBreaks(Philadelphia,"pop_Change"),1,5)) +
    scale_colour_manual(values = c("black","red")) + 
    labs(title = "Projected Population, 2020: Montgomery") + mapTheme() +
    guides(fill = guide_legend(order = 1), colour = guide_legend(order = 2)), ncol=2)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

# Chester
Philadelphia <-
  dat2 %>%
  mutate(Development_Demand = predict(Model6, dat2, type="response")) %>%
  filter(NAME == "Chester") 

Philadelphia_landUse <- rbind(
  filter(Philadelphia, forest11 == 1 | wetlands11 == 1 ) %>%
    dplyr::select() %>% mutate(Land_Use = "Not Suitable"),
  filter(Philadelphia, developed11 == 1) %>%
    dplyr::select() %>% mutate(Land_Use = "Developed"))

grid.arrange(
  ggplot() +
    geom_sf(data=Philadelphia, aes(fill=factor(ntile(Development_Demand,5))), colour=NA) +
    geom_point(data=Philadelphia_landUse, aes(x=xyC(Philadelphia_landUse)[,1], 
                                              y=xyC(Philadelphia_landUse)[,2], colour=Land_Use),
               shape = 15, size = 2) +
    geom_sf(data=st_intersection(DVRHighways,filter(studyAreaCounties, NAME=="Chester")), size=2) +
    scale_fill_manual(values = palette5, name="Development\nDemand",
                      labels=substr(quintileBreaks(Philadelphia,"Development_Demand"),1,5)) +
    scale_colour_manual(values = c("black","red")) + 
    labs(title = "Development Potential, 2020: Chester") + mapTheme() +
    guides(fill = guide_legend(order = 1), colour = guide_legend(order = 2)),
  
  ggplot() +
    geom_sf(data=Philadelphia, aes(fill=factor(ntile(pop_Change,5))), colour=NA) +
    geom_point(data=Philadelphia_landUse, aes(x=xyC(Philadelphia_landUse)[,1], 
                                              y=xyC(Philadelphia_landUse)[,2], colour=Land_Use),
               shape = 15, size = 2) +
    geom_sf(data=st_intersection(DVRHighways,filter(studyAreaCounties, NAME=="Chester")), size=2) +
    scale_fill_manual(values = palette5, name="Population\nChange",
                      labels=substr(quintileBreaks(Philadelphia,"pop_Change"),1,5)) +
    scale_colour_manual(values = c("black","red")) + 
    labs(title = "Projected Population, 2020: Chester") + mapTheme() +
    guides(fill = guide_legend(order = 1), colour = guide_legend(order = 2)), ncol=2)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

# Mercer

Philadelphia <-
  dat2 %>%
  mutate(Development_Demand = predict(Model6, dat2, type="response")) %>%
  filter(NAME == "Mercer") 

Philadelphia_landUse <- rbind(
  filter(Philadelphia, forest11 == 1 | wetlands11 == 1 ) %>%
    dplyr::select() %>% mutate(Land_Use = "Not Suitable"),
  filter(Philadelphia, developed11 == 1) %>%
    dplyr::select() %>% mutate(Land_Use = "Developed"))

grid.arrange(
  ggplot() +
    geom_sf(data=Philadelphia, aes(fill=factor(ntile(Development_Demand,5))), colour=NA) +
    geom_point(data=Philadelphia_landUse, aes(x=xyC(Philadelphia_landUse)[,1], 
                                              y=xyC(Philadelphia_landUse)[,2], colour=Land_Use),
               shape = 15, size = 2) +
    geom_sf(data=st_intersection(DVRHighways,filter(studyAreaCounties, NAME=="Mercer")), size=2) +
    scale_fill_manual(values = palette5, name="Development\nDemand",
                      labels=substr(quintileBreaks(Philadelphia,"Development_Demand"),1,5)) +
    scale_colour_manual(values = c("black","red")) + 
    labs(title = "Development Potential, 2020: Mercer") + mapTheme() +
    guides(fill = guide_legend(order = 1), colour = guide_legend(order = 2)),
  
  ggplot() +
    geom_sf(data=Philadelphia, aes(fill=factor(ntile(pop_Change,5))), colour=NA) +
    geom_point(data=Philadelphia_landUse, aes(x=xyC(Philadelphia_landUse)[,1], 
                                              y=xyC(Philadelphia_landUse)[,2], colour=Land_Use),
               shape = 15, size = 2) +
    geom_sf(data=st_intersection(DVRHighways,filter(studyAreaCounties, NAME=="Mercer")), size=2) +
    scale_fill_manual(values = palette5, name="Population\nChange",
                      labels=substr(quintileBreaks(Philadelphia,"pop_Change"),1,5)) +
    scale_colour_manual(values = c("black","red")) + 
    labs(title = "Projected Population, 2020: Mercer") + mapTheme() +
    guides(fill = guide_legend(order = 1), colour = guide_legend(order = 2)), ncol=2)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

# Delaware

Philadelphia <-
  dat2 %>%
  mutate(Development_Demand = predict(Model6, dat2, type="response")) %>%
  filter(NAME == "Delaware") 

Philadelphia_landUse <- rbind(
  filter(Philadelphia, forest11 == 1 | wetlands11 == 1 ) %>%
    dplyr::select() %>% mutate(Land_Use = "Not Suitable"),
  filter(Philadelphia, developed11 == 1) %>%
    dplyr::select() %>% mutate(Land_Use = "Developed"))

grid.arrange(
  ggplot() +
    geom_sf(data=Philadelphia, aes(fill=factor(ntile(Development_Demand,5))), colour=NA) +
    geom_point(data=Philadelphia_landUse, aes(x=xyC(Philadelphia_landUse)[,1], 
                                              y=xyC(Philadelphia_landUse)[,2], colour=Land_Use),
               shape = 15, size = 2) +
    geom_sf(data=st_intersection(DVRHighways,filter(studyAreaCounties, NAME=="Delaware")), size=2) +
    scale_fill_manual(values = palette5, name="Development\nDemand",
                      labels=substr(quintileBreaks(Philadelphia,"Development_Demand"),1,5)) +
    scale_colour_manual(values = c("black","red")) + 
    labs(title = "Development Potential, 2020: Delaware") + mapTheme() +
    guides(fill = guide_legend(order = 1), colour = guide_legend(order = 2)),
  
  ggplot() +
    geom_sf(data=Philadelphia, aes(fill=factor(ntile(pop_Change,5))), colour=NA) +
    geom_point(data=Philadelphia_landUse, aes(x=xyC(Philadelphia_landUse)[,1], 
                                              y=xyC(Philadelphia_landUse)[,2], colour=Land_Use),
               shape = 15, size = 2) +
    geom_sf(data=st_intersection(DVRHighways,filter(studyAreaCounties, NAME=="Delaware")), size=2) +
    scale_fill_manual(values = palette5, name="Population\nChange",
                      labels=substr(quintileBreaks(Philadelphia,"pop_Change"),1,5)) +
    scale_colour_manual(values = c("black","red")) + 
    labs(title = "Projected Population, 2020: Delaware") + mapTheme() +
    guides(fill = guide_legend(order = 1), colour = guide_legend(order = 2)), ncol=2)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

# Camden
Philadelphia <-
  dat2 %>%
  mutate(Development_Demand = predict(Model6, dat2, type="response")) %>%
  filter(NAME == "Camden") 

Philadelphia_landUse <- rbind(
  filter(Philadelphia, forest11 == 1 | wetlands11 == 1 ) %>%
    dplyr::select() %>% mutate(Land_Use = "Not Suitable"),
  filter(Philadelphia, developed11 == 1) %>%
    dplyr::select() %>% mutate(Land_Use = "Developed"))

grid.arrange(
  ggplot() +
    geom_sf(data=Philadelphia, aes(fill=factor(ntile(Development_Demand,5))), colour=NA) +
    geom_point(data=Philadelphia_landUse, aes(x=xyC(Philadelphia_landUse)[,1], 
                                              y=xyC(Philadelphia_landUse)[,2], colour=Land_Use),
               shape = 15, size = 2) +
    geom_sf(data=st_intersection(DVRHighways,filter(studyAreaCounties, NAME=="Camden")), size=2) +
    scale_fill_manual(values = palette5, name="Development\nDemand",
                      labels=substr(quintileBreaks(Philadelphia,"Development_Demand"),1,5)) +
    scale_colour_manual(values = c("black","red")) + 
    labs(title = "Development Potential, 2020: Camden") + mapTheme() +
    guides(fill = guide_legend(order = 1), colour = guide_legend(order = 2)),
  
  ggplot() +
    geom_sf(data=Philadelphia, aes(fill=factor(ntile(pop_Change,5))), colour=NA) +
    geom_point(data=Philadelphia_landUse, aes(x=xyC(Philadelphia_landUse)[,1], 
                                              y=xyC(Philadelphia_landUse)[,2], colour=Land_Use),
               shape = 15, size = 2) +
    geom_sf(data=st_intersection(DVRHighways,filter(studyAreaCounties, NAME=="Camden")), size=2) +
    scale_fill_manual(values = palette5, name="Population\nChange",
                      labels=substr(quintileBreaks(Philadelphia,"pop_Change"),1,5)) +
    scale_colour_manual(values = c("black","red")) + 
    labs(title = "Projected Population, 2020: Camden") + mapTheme() +
    guides(fill = guide_legend(order = 1), colour = guide_legend(order = 2)), ncol=2)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

# Burlington

Philadelphia <-
  dat2 %>%
  mutate(Development_Demand = predict(Model6, dat2, type="response")) %>%
  filter(NAME == "Burlington") 

Philadelphia_landUse <- rbind(
  filter(Philadelphia, forest11 == 1 | wetlands11 == 1 ) %>%
    dplyr::select() %>% mutate(Land_Use = "Not Suitable"),
  filter(Philadelphia, developed11 == 1) %>%
    dplyr::select() %>% mutate(Land_Use = "Developed"))

grid.arrange(
  ggplot() +
    geom_sf(data=Philadelphia, aes(fill=factor(ntile(Development_Demand,5))), colour=NA) +
    geom_point(data=Philadelphia_landUse, aes(x=xyC(Philadelphia_landUse)[,1], 
                                              y=xyC(Philadelphia_landUse)[,2], colour=Land_Use),
               shape = 15, size = 2) +
    geom_sf(data=st_intersection(DVRHighways,filter(studyAreaCounties, NAME=="Burlington")), size=2) +
    scale_fill_manual(values = palette5, name="Development\nDemand",
                      labels=substr(quintileBreaks(Philadelphia,"Development_Demand"),1,5)) +
    scale_colour_manual(values = c("black","red")) + 
    labs(title = "Development Potential, 2020: Burlington") + mapTheme() +
    guides(fill = guide_legend(order = 1), colour = guide_legend(order = 2)),
  
  ggplot() +
    geom_sf(data=Philadelphia, aes(fill=factor(ntile(pop_Change,5))), colour=NA) +
    geom_point(data=Philadelphia_landUse, aes(x=xyC(Philadelphia_landUse)[,1], 
                                              y=xyC(Philadelphia_landUse)[,2], colour=Land_Use),
               shape = 15, size = 2) +
    geom_sf(data=st_intersection(DVRHighways,filter(studyAreaCounties, NAME=="Burlington")), size=2) +
    scale_fill_manual(values = palette5, name="Population\nChange",
                      labels=substr(quintileBreaks(Philadelphia,"pop_Change"),1,5)) +
    scale_colour_manual(values = c("black","red")) + 
    labs(title = "Projected Population, 2020: Burlington") + mapTheme() +
    guides(fill = guide_legend(order = 1), colour = guide_legend(order = 2)), ncol=2)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

# Gloucester

Philadelphia <-
  dat2 %>%
  mutate(Development_Demand = predict(Model6, dat2, type="response")) %>%
  filter(NAME == "Gloucester") 

Philadelphia_landUse <- rbind(
  filter(Philadelphia, forest11 == 1 | wetlands11 == 1 ) %>%
    dplyr::select() %>% mutate(Land_Use = "Not Suitable"),
  filter(Philadelphia, developed11 == 1) %>%
    dplyr::select() %>% mutate(Land_Use = "Developed"))

grid.arrange(
  ggplot() +
    geom_sf(data=Philadelphia, aes(fill=factor(ntile(Development_Demand,5))), colour=NA) +
    geom_point(data=Philadelphia_landUse, aes(x=xyC(Philadelphia_landUse)[,1], 
                                              y=xyC(Philadelphia_landUse)[,2], colour=Land_Use),
               shape = 15, size = 2) +
    geom_sf(data=st_intersection(DVRHighways,filter(studyAreaCounties, NAME=="Gloucester")), size=2) +
    scale_fill_manual(values = palette5, name="Development\nDemand",
                      labels=substr(quintileBreaks(Philadelphia,"Development_Demand"),1,5)) +
    scale_colour_manual(values = c("black","red")) + 
    labs(title = "Development Potential, 2020: Gloucester") + mapTheme() +
    guides(fill = guide_legend(order = 1), colour = guide_legend(order = 2)),
  
  ggplot() +
    geom_sf(data=Philadelphia, aes(fill=factor(ntile(pop_Change,5))), colour=NA) +
    geom_point(data=Philadelphia_landUse, aes(x=xyC(Philadelphia_landUse)[,1], 
                                              y=xyC(Philadelphia_landUse)[,2], colour=Land_Use),
               shape = 15, size = 2) +
    geom_sf(data=st_intersection(DVRHighways,filter(studyAreaCounties, NAME=="Gloucester")), size=2) +
    scale_fill_manual(values = palette5, name="Population\nChange",
                      labels=substr(quintileBreaks(Philadelphia,"pop_Change"),1,5)) +
    scale_colour_manual(values = c("black","red")) + 
    labs(title = "Projected Population, 2020: Gloucester") + mapTheme() +
    guides(fill = guide_legend(order = 1), colour = guide_legend(order = 2)), ncol=2)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

  1. Supply-side Prediction with new highway
NEWHighways <-
  read_sf("new_highways.shp") %>%
  st_transform(st_crs(DVR)) %>%
  st_intersection(DVR)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
plot(NEWHighways)

ggplot() +
  geom_point(data=fishnet, 
             aes(x=xyC(fishnet)[,1], y=xyC(fishnet)[,2],colour=lc_change),size=1.5) +
  geom_sf(data=NEWHighways) +
  scale_colour_manual(values = palette2,na.translate = F,
                      labels=c("No Change","New Development")) +
  labs(title = "New Development and Highways",
       subtitle = "As fishnet centroids") +
  mapTheme()
## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x
## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

emptyRaster2 <- lc_change
emptyRaster2[] <- NA


NEWhighway_raster <- 
  as(NEWHighways,'Spatial') %>%
  rasterize(.,emptyRaster)

NEWhighway_raster_distance <- distance(NEWhighway_raster)
names(NEWhighway_raster_distance) <- "new_distance_highways"

NEWhighwayPoints <-
  rasterToPoints(NEWhighway_raster_distance) %>%
  as.data.frame() %>%
  st_as_sf(coords = c("x", "y"), crs = st_crs(DVR_fishnet))


NEWhighwayPoints_fishnet <- 
  aggregate(NEWhighwayPoints, DVR_fishnet, mean) %>%
  mutate(new_distance_highways = ifelse(is.na(new_distance_highways),0,new_distance_highways))

ggplot() +
  geom_sf(data=DVR) +
  geom_point(data=NEWhighwayPoints_fishnet, aes(x=xyC(NEWhighwayPoints_fishnet)[,1], 
                                             y=xyC(NEWhighwayPoints_fishnet)[,2], 
                                             colour=factor(ntile(new_distance_highways,5))),size=1.5) +
  scale_colour_manual(values = palette5,
                      labels=substr(quintileBreaks(NEWhighwayPoints_fishnet,"new_distance_highways"),1,8),
                      name="Quintile\nBreaks") +
  geom_sf(data=NEWHighways, colour = "red") +
  labs(title = "Distance to NEW Highways",
       subtitle = "As fishnet centroids; Highways visualized in red") +
  mapTheme()
## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

## Warning in st_centroid.sf(aPolygonSF): st_centroid assumes attributes are
## constant over geometries of x

dat2 <- cbind( fishnet, NEWhighwayPoints_fishnet, fishnetPopulation, aggregatedRasters) %>% dplyr::select(lc_change, developed, forest, farm, wetlands, otherUndeveloped, water, pop_2000, pop_2010, pop_Change, new_distance_highways,lagDevelopment) %>% st_join(studyAreaCounties) %>% mutate(developed10 = ifelse(lc_change == 1 & developed == 1, 0, developed)) %>% filter(water == 0)

dat2 %>% mutate(new_distance_highways = replace_na(new_distance_highways, 0))

dat2 %>% mutate(lc_change = replace_na(lc_change, 0))

dat2 %>% dplyr::select(new_distance_highways,lagDevelopment,lc_change) %>% gather(Variable, Value, -lc_change, -geometry) %>% ggplot(., aes(lc_change, Value, fill=lc_change,), Variable[!is.na(Variable$lc_change),]) + geom_bar(position = “dodge”, stat = “summary”, fun.y = “mean”) + facet_wrap(~Variable) + scale_fill_manual(values = palette2,na.translate = F, labels=c(“No Change”,“New Development”), name="“) + labs(title=”New development as a function of the countinuous variables") + plotTheme()

dat2 %>% dplyr::select(pop_2000,pop_2010,pop_Change,lc_change) %>% gather(Variable, Value, -lc_change, -geometry) %>% ggplot(., aes(lc_change, Value, fill=lc_change)) + geom_bar(position = “dodge”, stat = “summary”, fun.y = “mean”) + facet_wrap(~Variable) + scale_fill_manual(values = palette2,na.translate = F, labels=c(“No Change”,“New Development”), name="“) + labs(title=”New development as a function of factor variables") + plotTheme()

dat2 %>% dplyr::select(lc_change:otherUndeveloped,developed) %>% gather(Land_Cover_Type, Value, -lc_change, -geometry) %>% st_set_geometry(NULL) %>% group_by(lc_change, Land_Cover_Type) %>% summarize(n = sum(as.numeric(Value))) %>% ungroup() %>% mutate(Conversion_Rate = paste0(round(100 * n/sum(n), 2), “%”)) %>% filter(lc_change == 1) %>% dplyr::select(Land_Cover_Type,Conversion_Rate) %>% kable() %>% kable_styling(full_width = F)

set.seed(3456) trainIndex2 <- createDataPartition(dat2$developed, p = .50, list = FALSE, times = 1) datTrain2 <- dat2[ trainIndex2,] datTest2 <- dat2[-trainIndex2,]

nrow(dat2)

Model11 <- glm(lc_change ~ wetlands + forest + farm + otherUndeveloped, family=“binomial”(link=“logit”), data = datTrain2)

Model12 <- glm(lc_change ~ wetlands + forest + farm + otherUndeveloped + lagDevelopment, family=“binomial”(link=“logit”), data = datTrain2)

Model13 <- glm(lc_change ~ wetlands + forest + farm + otherUndeveloped + lagDevelopment + pop_2000, family=“binomial”(link=“logit”), data = datTrain2)

Model14 <- glm(lc_change ~ wetlands + forest + farm + otherUndeveloped + lagDevelopment + pop_2000 + pop_2010, family=“binomial”(link=“logit”), data = datTrain2)

Model15 <- glm(lc_change ~ wetlands + forest + farm + otherUndeveloped + lagDevelopment + pop_Change, family=“binomial”(link=“logit”), data = datTrain2)

Model16 <- glm(lc_change ~ wetlands + forest + farm + otherUndeveloped + lagDevelopment + pop_Change + new_distance_highways, family=“binomial”(link=“logit”), data = datTrain2)

modelList <- paste0(“Model”, 11:16) map_dfc(modelList, function(x)pR2(get(x)))[4,] %>% setNames(paste0(“Model”,11:16)) %>% gather(Model,McFadden) %>% ggplot(aes(Model,McFadden)) + geom_bar(stat=“identity”) + labs(title= “McFadden R-Squared by Model”) + plotTheme()

testSetProbs2 <- data.frame(class2 = datTest2$lc_change, probs2 = predict(Model16, datTest2, type=“response”))

ggplot(testSetProbs2, aes(probs2)) + geom_density(aes(fill=class2), alpha=0.5) + scale_fill_manual(values = palette2,na.translate = F, labels=c(“No Change”,“New Development”)) + labs(title = “Histogram of test set predicted probabilities”, x=“Predicted Probabilities”,y=“Density”) + plotTheme()

options(yardstick.event_first = FALSE)

testSetProbs2 <- testSetProbs2 %>% mutate(predClass2_05 = as.factor(ifelse(testSetProbs2\(probs2 >= 0.05 ,1,0)), predClass2_17 = as.factor(ifelse(testSetProbs2\)probs2 >= 0.17 ,1,0)))

testSetProbs2 %>% dplyr::select(-probs2) %>% gather(Variable, Value, -class2) %>% group_by(Variable) %>% summarize(Sensitivity = round(yardstick::sens_vec(class2,factor(Value)),2), Specificity = round(yardstick::spec_vec(class2,factor(Value)),2), Accuracy = round(yardstick::accuracy_vec(class2,factor(Value)),2)) %>% kable() %>% kable_styling(full_width = F)

predsForMap2 <-
dat2 %>% mutate(probs2 = predict(Model16, dat2, type=“response”) , Threshold_5_Pct = as.factor(ifelse(probs2 >= 0.05 ,1,0)), Threshold_17_Pct = as.factor(ifelse(probs2 >= 0.17 ,1,0))) %>% dplyr::select(lc_change,Threshold_5_Pct,Threshold_17_Pct) %>% gather(Variable,Value, -geometry) %>% st_cast(“POLYGON”)

ggplot() + geom_point(data=predsForMap2, aes(x=xyC(predsForMap2)[,1], y=xyC(predsForMap2)[,2], colour=Value)) + facet_wrap(~Variable) + scale_colour_manual(values = palette2, labels=c(“No Change”,“New Development”),na.translate = F, name="“) + labs(title=”New Development predictions - Low threshold") + mapTheme()

ConfusionMatrix.metrics2 <- dat2 %>% mutate(probs2 = predict(Model16, dat2, type=“response”) , Threshold_5_Pct = as.factor(ifelse(probs2 >= 0.05 ,1,0)), Threshold_17_Pct = as.factor(ifelse(probs2 >= 0.17 ,1,0))) %>% mutate(TrueP_05 = ifelse(lc_change == 1 & Threshold_5_Pct == 1, 1,0), TrueN_05 = ifelse(lc_change == 0 & Threshold_5_Pct == 0, 1,0), TrueP_17 = ifelse(lc_change == 1 & Threshold_17_Pct == 1, 1,0), TrueN_17 = ifelse(lc_change == 0 & Threshold_17_Pct == 0, 1,0)) %>% dplyr::select(., starts_with(“True”)) %>% gather(Variable, Value, -geometry) %>% st_cast(“POLYGON”)

ggplot(data=ConfusionMatrix.metrics2) + geom_point(aes(x=xyC(ConfusionMatrix.metrics2)[,1], y=xyC(ConfusionMatrix.metrics2)[,2], colour = as.factor(Value))) + facet_wrap(~Variable) + scale_colour_manual(values = palette2, labels=c(“Correct”,“Incorrect”),na.translate = F, name="“) + labs(title=”Development predictions - Low threshold") + mapTheme()

spatialCV <- function(dataFrame, uniqueID, dependentVariable, modelName) {

#initialize a data frame endList <- list()

#create a list that is all the spatial group unqiue ids in the data frame (ie counties)
uniqueID_List <- unique(dataFrame[[uniqueID]])
x <- 1 y <- length(uniqueID_List)

#create a counter and while it is less than the number of counties…
while(x <= y) { #call a current county
currentUniqueID <- uniqueID_List[x] #create a training set comprised of units not in that county and a test set of units #that are that county training2 <- dataFrame[ which(dataFrame[[uniqueID]] != uniqueID_List[x]),] testing2 <- dataFrame[ which(dataFrame[[uniqueID]] == uniqueID_List[x]),] #create seperate xy vectors trainingX2 <- training2[ , -which(names(training) %in% c(dependentVariable))] testingX2 <- testing2[ , -which(names(testing) %in% c(dependentVariable))]

trainY2 <- training2[[dependentVariable]] testY2 <- testing2[[dependentVariable]] #Calculate predictions on the test county as part of a data frame including the observed #outcome and the unique county ID
thisPrediction2 <- data.frame(class = testY, probs = predict(modelName, testingX2, type=“response”), county = currentUniqueID)

#Row bind the predictions to a data farme endList2 <- rbind(endList2, thisPrediction2) #iterate counter
x <- x + 1 } #return the final list of counties and associated predictions
return (as.data.frame(endList)) }

spatialCV_counties2 <- spatialCV(dat2,“NAME”,“lc_change”, Model16) %>% mutate(predClass2 = as.factor(ifelse(probs2 >= 0.17 ,1,0)))

spatialCV_metrics2 <- spatialCV_counties2 %>% group_by(county) %>% summarize(Observed_Change = sum(as.numeric(as.character(class))), Sensitivity = round(yardstick::sens_vec(class,predClass2),2), Specificity = round(yardstick::spec_vec(class,predClass2),2), Accuracy = round(yardstick::accuracy_vec(class,predClass2),2))

spatialCV_metrics2 %>% kable() %>% kable_styling(full_width = F)

dat2 <- dat2 %>% mutate(lagDevelopment = nn_function(xyC(.), xyC(filter(.,developed10 == 1)),2))

countyPopulation_2020 <- data.frame( NAME = c(“Bucks”,“Montgomery”,“Chester”,“Philadelphia”,“Delaware”, “Mercer”,“Camden”,“Gloucester”,“Burlington”), county_projection_2020 = c(640495, 840934, 543702, 1594787, 568337, 377328 , 514006, 307766, 459344)) %>% left_join( dat %>% st_set_geometry(NULL) %>% group_by(NAME) %>% summarize(county_population_2010 = round(sum(pop_2010))))

countyPopulation_2020 %>% gather(Variable,Value, -NAME) %>% ggplot(aes(reorder(NAME,-Value),Value)) + geom_bar(aes(fill=Variable), stat = “identity”) + scale_fill_manual(values = palette2, labels=c(“2020”,“2010”), name=“Population”) + labs(title=“Population Change by county: 2010 - 2020”, x=“County”, y=“Population”) + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + plotTheme()

dat_infill2 <- dat2 %>% #calculate population change left_join(countyPopulation_2020) %>% mutate(proportion_of_county_pop = pop_2010 / county_population_2010, pop_2020.infill = proportion_of_county_pop * county_projection_2020, pop_Change = round(pop_2020.infill - pop_2010),2) %>% dplyr::select(-county_projection_2020, -county_population_2010, -proportion_of_county_pop, -pop_2020.infill) %>% #predict for 2020 mutate(predict_2020.infill = predict(Model16,. , type=“response”)) #model6 not model16 dat_infill2 %>% ggplot() +
geom_point(aes(x=xyC(dat_infill2)[,1], y=xyC(dat_infill2)[,2], colour = factor(ntile(predict_2020.infill,5)))) + scale_colour_manual(values = palette5, labels=substr(quintileBreaks(dat_infill2,“predict_2020.infill”),1,4), name=“Quintile”) + geom_sf(data=studyAreaCounties, fill=NA, colour=“black”, size=1.5) + labs(title= “Development Demand in 2020: Predicted Probabilities”) + mapTheme()